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FOREWORD 

The  Operations  Research  Center  at  the  Massachusetts  Institute  of 
Technology  is  an  interdepartmental  activity  devoted  to  graduate  educa- 
tion and  research  in  the  field  of  operations  research.  The  work  of  the 
Center  is  supported,  in  part,  by  government  grants  and  contracts.  The 
work  reported  herein  was  supported  (in  part)  by  the  Office  of  Naval  Re- 
search under  Contract  N00014-75-C-0556  and  by  the  Department  of  Trans- 
portation Advanced  Research  Program  under  contract  DOT-TSC-1058. 

Richard  C.  Larson 

Jeremy  F.  Shapiro 

Co-Directors 


ABSTRACT 

Transportation  planning  plays  an  essential  role  in  shaping  regional  and 
urban  lifestyle.  Complex  decisions  regarding  policy  alternatives  for  rail- 
roads, shipping,  airline,  and  roadway  traffic  can  often  be,  and  often  have 
been,  analyzed  using  network  optimization  techniques.  In  this  paper,  we  survey 
applications  of  network  algorithms  to  transportation  planning,  stressing  net- 
work models  and  their  efficient  computer  implementation.  We  discuss  recent 
contributions  concerning  shortest  paths,  minimum  cost  network  flows,  traffic 
equilibrium,  vehicle  routing,  and  network  design  and  we  enumerate  several  open 
research  problems.  Much  of  our  discussion  reflects  an  emerging  theme  in  the 
analysis  of  transportation  problems,  the  blending  of  Ideas  from  transportation 
science,  computer  science,  and  operations  research. 
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1.  Introduction;  Modeling  and  Implementation  Issues 

Transportation  is  vital  to  any  society.  It  exerts  great  Influence  on 
the  flow  of  goods,  services,  and  information,  on  the  location  of  homes  and 
industry,  on  the  use  of  recreational  and  cultural  activities;  in  many  ways,  it 
does  much  to  define  the  character  of  our  lives. 

Needless  to  say,  planning  for  effective  transportation  is  a complicated 
process  requiring  estimation  of  transportation  needs,  technological  Innovations, 
assessment  of  new  and  proposed  investments,  and  efficient  management  of  exist- 
ing facilities.  In  most  planning  efforts,  it  is  natural  to  view  a transportation 
system  as  a transportation  network  with  a number  of  nodes  (representing  street 
intersections,  depots,  ports,  cities,  and  so  on)  and  a number  of  links,  arcs,  or 
edges  (e.g.,  streets,  air  and  ship  routes,  or  subway  channels).  Network  models, 
hence,  become  the  focal  point  for  a great  deal  of  analysis  of  transportation 
systems.  The  models  may  be  normative,  the  optimization  of  system  performance 
(usually  with  costs  or  travel  times  associated  with  the  links)  being  the 
objective.  Or,  they  may  be  predictive.  How  will  users  (individuals,  private 
and  public  organizations)  and  the  transportation  industry  Itself  respond  to 
various  policy  alternatives?  We  consider  both  types  of  models  in  this  paper. 

Much  of  the  recent  research  in  netwrk  analysis  for  transportation 
planning  has  involved  a blending  of  ideas  from  transportation  science,  opera- 
tions research,  and  computer  science.  The  Interplay  between  modeling,  algo- 
rithms, and  their  efficient  computer  implementation  has  become  a dominant  theme. 
To  illustrate  the  nature  of  this  research,  and  at  the  same  time  consider  areas 
rich  in  applications,  we  have  selected  the  following  topics  for  discussion: 
shortest  paths,  minimum  cost  network  flows,  traffic  equilibrium,  vehicle  rout- 
ing, and  network  design. 

We  might  note  that  our  coverage  of  the  first  three  of  these  topics, 
when  contrasted  with  Potts  and  Oliver's  highly-regarded  book [105]  published 
only  a few  years  ago,  is  Indicative  of  the  current  level  of  activity  and 
recent  progress  in  transportation  planning  and  network  analysis. 

Transportation  Models 

The  streets  of  a city  form  an  obvious  network  where  nodes  (numbered  for 
identification  purposes)  represent  locations  and  Intersections  on  roads,  and 
links  represent  the  roads  themselves.  We  distinguish  between  one-way  and  two- 
way  streets  by  using  a directed  link  for  the  former  and  an  undirected  link  (or 
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two  directed  links)  for  Che  laCter.  As  an  Illustration  of  how  a real  trans- 
portation system  may  be  abstracted  by  a network  model.  Figure  1 shows  a segment 
of  a fictitious  city  street  network. 


In  dealing  with  a network  model  of  a real  transportation  system,  trans- 
1 porCation  planners  typically  associate  various  parameters  with  the  nodes  and 

links.  For  example,  each  link  of  an  urban  road  network  may  have  values  for  the 

I 

following  Items: 

(1)  number  of  traffic  lanes, 

(11)  road  length, 

(111)  average  travel  time, 

(Iv)  average  vehicle  speeds, 

(v)  average  dally  traffic  flow, 

(vl)  peak  hour  flows, 

(vll)  capacity 

(vlll)  total  monetary  cost  (Including  tolls) . 

These  values  are  frequently  combined  In  order  to  obtain  a single  measure  of  cost 
or  distance  on  the  link. 

Different  types  of  networks  arise  In  other  settings.  What,  for  example, 
la  the  maxlouffl  Income  for  an  airline  system?  Given  a number  of  possible  non- 
stop services  or  flights,  an  associated  expected  Income  for  each  service,  and 
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an  overnight  holding  cost  for  an  aircraft,  what  set  of  services  should  be 
flown  for  maximum  system  Income  (Simpson  [ll^ describes  this  model  and  many 
others).  Figure  2 provides  a representation  of  such  a problem  as  a time-space 
network.  There  are  three  geographical  locations  A,  B,  and  C.  The  network  con- 
sists of  nodes  which  Indicate  both  geographic  location  and  time  of  day.  Poten- 
tial flights  are  shown  by  "service”  arcs  joining  geographic  locations  at 
various  times  of  the  day;  expected  Incomes  are  associated  with  these  arcs. 

An  arc  from  the  end  of  the  day  to  the  beginning  of  the  next  day  corresponds  to 
holding  an  aircraft  overnight.  There  are  dally  rental  costs  associated  with 
these  arcs.  The  problem  Is  to  find  the  maximum  revenue  route  schedule  subject 
to  the  capacity  limitation  that  at  most  one  plane  flies  any  service  arc. 
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Although  one  might  usually  associate  vehicles  with  the  examples  of 
Figures  1 and  2,  networks  model  other  aspects  of  transportation  planning  as 
well,  such  as  the  flow  of  passengers,  cargo,  and  vehicle  crews.  In  fact,  one 
of  the  most  noteworthy  features  of  transportation  systems,  and  their  representa- 
tion as  models.  Is  that  they  usually  Involve  several  different  commodities.  In 
many  Instances  commodities  will  be  distinguished  by  their  points  of  origin  and 
destination.  Passengers  traveling  from  New  York  to  Los  Angeles  are  not  indis- 
tinguishable from  those  traveling  from  Washington  to  San  Francisco,  even  though 
they  may  share  some  of  the  same  transportation  facilities;  otherwise,  the  model 
could  route  the  New  York  passengers  to  San  Francisco  and  those  from  Washington 
to  Los  Angeles.  When  passengers  (or  cargo  or  crews)  constitute  one  type  of 
commodity  and  vehicles  another,  the  models  are  further  complicated  because  the 
vehicle  flows  define  possible  routes  and  capacity  limits  for  passenger  travel. 

In  any  event,  realistic  modeling  of  transportation  systems  often  results  In 
multicommodity  models. 

The  taxonomy  of  strategic,  tactical,  and  operational  decision  making, 
as  outlined  In  Table  1,  helps  to  distinguish  between  different  types  of  models 
for  transportation  planning.  In  our  more  detailed  discussion  of  particular 
models,  we  shall  consider  problems  from  each  general  category. 

Implementing  Transportation  Models 

Since  transportation  models,  like  most  others,  need  to  be  solved  repeat- 
edly In  order  to  study  modeling  assumptions,  to  perform  sensitivity  analysis, 
and  to  address  changes  over  time.  It  becomes  essential  that  algorithms  be 
designed  to  run  efficiently.  For  networks,  the  manner  In  which  data  Is  stored 
and  manipulated  often  has  a significant  Impact  upon  an  algorithm's  performance. 
Also,  because  of  the  nature  of  network  topology,  special  techniques  are  available 
to  structure  problem  data  within  a computer. 

Suppose  that  we  need  to  store  a network  with  a constant  per  unit  cost 
for  each  arc.  Perhaps  the  easiest  scheme  to  work  with, and  yet  the  most  Inef- 
ficient In  terms  of  storage  conservation, Is  a matrix  representation  which  has 
as  the  l,j  th  entry  the  cost  of  the  arc  from  1 to  j,  or  • If  no  arc  exists.  If 
the  network  has  n nodes  and  E arcs,  then  n^  locations  are  required.  Note  that 
for  sparse  networks  most  entries  will  be  ••.  Another  way  of  storing  network  data 

1 is  known  as  the  "ladder  representation."  For  each  arc  we  record  its  origin 

node,  its  destination  node,  and  its  cost.  This  approach  calls  for  3E  locations. 

i 
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Table  1.  Taxonomy  of  Transportation  Planning  Models 
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A third  representation,  the  "forward  star  representation,"  records  the  arcs 
ordered  by  origin  node.  An  arc  list  contains  for  arc  k,  its  destination  node 
and  its  cost.  An  auxiliary  node  list  records  for  node  i,  the  first  entry  in 
the  arc  list  originating  from  that  node.  The  scheme  requires  n + 2E  storage 
locations. 

Now  suppose  the  cost  functions  are  of  a more  complex  nature.  For  example, 

2 

if  the  cost  on  an  arc  with  flow  x is  given  by  a quadratic  function  ax  + bx  + c, 

then  we  could  either  keep  three  cost  matrices,  or  store  three  cost  vectors  A, 

B,  and  C using  the  ladder  or  forward  star  representations.  These  approaches 
2 

would  require  3n  , 5E,  and  n + 4E  storage  locations  respectively.  The  storage 
schemes  would  be  used  in  the  same  way  to  record  other  data,  such  as  arc  capacities. 

Depending  on  the  functional  form  of  the  cost  functions,  the  sparsity  of  the 
network,  and  the  amount  of  data  manipulation  required  by  an  algorithm,  the  user 
must  determine  the  best  network  representation  for  his  particular  application. 

We  discuss  this  issue,  and  how  it  relates  to  implementing  network  algorithms, 
in  subsequent  sections  of  this  paper. 

For  additional  general  information  regarding  transportation  networks,  the 
excellent  surveys  by  Bradley  [19]  and  Gazis  [56]  are  recommended.  Also,  see 
Gartner  et  al.  [53]  and  Steenbrink  [118].  For  an  extensive  bibliography  on 
network  optimization,  see  Golden  and  Magnanti  [68]. 
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2.  Shortest  Paths 

Despite  the  number  of  papers  on  shortest  path  problems  surveyed  by 
Dreyfus  [40]  and  later  by  Gllslnn  and  Wltzgall  [65],  new  insights  regarding 
this  class  of  problems  continue  to  emerge.  In  the  last  few  years,  a number  of 
algorithmic  Improvements  have  been  reported  which  impact  directly  on  transpor- 
tation planning.  Tn  this  section,  we  outline  some  of  these  recent  contribu- 
tions. 

Overview 

Shortest  path  problems  are  pervasive  in  transportation  planning  for 
several  reasons.  One  of  the  primary  objectives  of  any  traveler  (a  passenger 
or  a carrier)  is  to  move  from  one  point  a to  another  point  b,  along  a shortest, 
cheapest . or  most  comfortable  path.  Associating  flow  costs  (distances  or  com- 
fort factors)  with  arcs  in  a network,  the  traveler  seeks  the  minimum  cost  path 
from  a to  b.  In  economic  terms,  there  is  a supply  of  one  or  more  units  at  node 
a,  a demand  for  these  units  at  node  b,  and  a link  flow  cost  assigned  to  each 
link  in  the  network. 

Shortest  path  problems  also  arise  in  situations  where  this  model,  by 
Itself,  is  not  appropriate,  such  as  when  the  route  selected  by  one  traveler 
effects  the  cost  of  routes  taken  by  other  travelers.  In  this  paper,  we  discuss 
a number  of  important  problems  and  techniques  in  network  optimization  relating 
to  transportation.  In  some  way,  each  problem  relies  or  builds  upon  a shortest 
path  algorithm.  The  minimum  cost  network  flow  problem  is  a generalization  permitting 
supplies  and  demands  for  flow  at  various  points  in  the  network  and  flow  capaci- 
ties on  the  links.  The  network  design  problem  Introduces  link  construction  pos- 
sibilities. When  urban  transportation  planners  try  to  forecast  traffic,  the 
shortest  path  problem  becomes  an  Important  subproblem.  The  vehicle  routing 
problem  requires  a shortest  path  matrix  as  input.  As  these  problems  illustrate, 
a shortest  path  algorithm  is  at  the  core  of  many  problems  in  transportation 
planning. 

In  terms  of  modeling  transportation  networks,  it  is  Important  to  realize 
that,  in  computing  shortest  paths,  total  travel  time  between  points  a and  b 
depends  not  only  on  link  travel  times,  but  also  on  delays  at  intersections, 
often  attributable  to  left  hand  turns.  Network  formulations  model  these  situations 
by  imposing  turn  penalties,  that  is,  by  associating  costs  or  delays  with  turns 
at  nodes.  Turn  prohibitions,  which  are  enforced  as  policies  in  many  transporta- 
tion systems,  can  be  regarded  as  turns  with  infinite  penalties. 
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Several  researchers  have  proposed  algorithms  for  determining  shortest 
routes  in  networks  with  turn  penalties  (see  Potts  and  Oliver  U05]  and  Kirby  and 
Potts  [88]  for  details).  Although  we  do  not  pursue  this  topic  here,  we  mention 
this  issue  because  of  its  important  modeling  implications.  An  example  other 
than  a road  network  that  might  be  modeled  with  turn  penalties  is  a subway  system 
with  many  different  lines.  Switching  lines  involves  a delay  and  possibly  a 
transfer  charge. 

Transportation  planning  is  not  the  only  setting  in  which  shortest  path 
problems  are  of  interest.  Similar  applications  arise  in  computer-communication 
studies.  In  addition,  shortest  path  problems  often  become  subproblems  for  more 
complex  problems  such  as  in  group  theoretic  integer  prograianing  (see  Shapiro  [ll3], 
Chen  and  Zlonts  [24],  Frieze  [52],  and  Denardo  and  Fox  [34]).  In  fact,  computa- 
tional studies  of  shortest  path  algorithms  have  inspired  research  in  sorting, 
data  structures,  and  list  processing  by  operations  researchers  and  computer 
scientists  alike. 

For  a given  network  G “ (N,  A,  D)  with  node  set  N,  arc  set  A,  and  arc 
costs  given  by  the  matrix  D = [d(i ,j )],  there  are  five  shortest  path  problems  of 
general  interest. 

(1)  Find  the  shortest  path  from  a specific  origin  s to  a specific  destina- 
tion t; 

(2)  Find  the  shortest  paths  from  a specific  origin  s to  all  other  nodes; 

(3)  Find  the  shortest  paths  between  all  pairs  of  nodes; 

(4)  Find  the  shortest  path  between  an  origin-destination  pair  that  passes 

through  specified  nodes; 

(5)  Find  the  second,  third,  and  so  on,  shortest  paths. 

The  distance  entries  d(l,J)  can  be  positive,  negative,  or  zero  provided  that 
there  exists  no  cycle  whose  total  cost  is  negative.  If  a negative  cycle  did 
exist,  costs  would  be  minimized  by  traversing  it  infinitely  often. 

Implementation  Issues 

Because  shortest  path  problems  are  so  central  to  transportation  science, 
efficient  implementation  of  computer  codes  for  these  problems  often  translates 
into  substantial  savings.  At  times,  the  efficiency  of  a code  dictates  the  size 
of  networks,  and  hence  the  detail  of  modeling,  that  can  be  analyzed.  Furthermore, 

In  real-time  planning  situations,  fast  computer  codes  become  a necessity. 

In  the  following  discussion,  we  focus  primarily  on  the  second  problem  listed 
above  which  Is,  perhaps,  the  most  common.  We  view  this  problem  and  Bellman's  algo- 
rithm [12]  for  solving  it  as  a vehicle  for  illustrating  how  rather  minor  changes  in  a 

- — * 
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code's  implementation  can  lead  to  significant  reductions  in  computer  running 
time.  Both  Bellman's  algorithm  [12],  and  a modification  of  It  proposed  by  Pape 
^04]  that  we  will  considei;  are  classified  as  "label-correcting"  procedures  (see 
[65]),  in  the  sense  that  tentative  shortest  path  distances  assigned  to  the  nodes 
are  revised  until  true  shortest  path  distances  are  determined.  We  outline  each 
procedure  below: 

BELLMAN'S  ALGORITHM  (also  known  as  the  Ford-Bellman-Moore  Procedure  [^0]) 
Definitions 

l(v)  is  the  length  of  the  current  "shortest  pathf' f rom  node  s to  node  v. 

p(v)  is  the  predecessor  of  v In  the  current  "shortest  path"to  this  node. 

d(l,k)  Is  Che  length  of  arc  (l,k)  e A. 

Initialization 

i(v)  . {0  if  V - s 
^ ^ ® otherwise 

p(v)  “ 0 for  all  V. 

node  s is  the  first  element  on  list  T. 

Basic  Computation 

Select  the  top  element  1 from  list  T.  For  every  node  k such  that 
(i,k}  c A,  perform  the  following  test: 

If  1(1)  + d(i,k)  < Kk)  then 

(a)  l(k)  - 1(1)  + d(i.k) 
p(k)  ■ 1,  and 

(b)  place  k at  bottom  of  T,  if  it  is  not  already  on  the  list. 
Reducing  the  List  Size 

Remove  (or  cross  out)  node  1 from  Che  list.  Terminate  Che  procedure 
if  the  list  T is  now  empty.  Otherwise  return  to  the  basic  computation. 

PAPE’S  ALGORITHM 

This  procedure  Is  the  same  as  the  previous  one  except  that  we  replace  (b) 
of  the  basic  computation  step  with: 

(b')  If  k Is  already  on  list  T,  do  not  add  It  again. 

If  k has  not  yet  been  on  the  list,  place  it  at  the  bottom  of  T. 

If  k has  already  been  processed  (that  is,  was  on  the  list  once  before 
but  Is  not  currently),  then  enter  k at  the  top  of  the  list. 

We  point  out  that  both  algorithms  are  based  on  the  following  fundamental 
recursion 

1(1)  min  {)l(J)  + d(J,l)} 

J 
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* where  the  labels  l(v)  are  updated  whenever  a path  of  one  additional  arc  has  a 

smaller  length  than  the  previous  best  path.  In  addition,  ench  algorithm 
requires  on  the  order  of  n^  additions  and  comparisons  in  the  worst  case,  where 
n is  the  number  of  nodes  in  the  network.  There  are  other  shortest  path  algo- 
rithms known  as  "label-setting"  procedures  ( see  [65])  that  only  require  on  the 
order  of  operations  in  the  worst  case. 

The  following  example  will  Illustrate  the  computational  advantage  of 
the  second  approach  over  the  first.  Furthermore,  recent  computational  studies 
by  Pape  [l04j  and  Klingman  et  al.  [89]  demonstrate  that  this  approach  seems  to 
outperform  other  types  of  shortest  path  algorithms  as  well.  Including  frequently 
advocated  "label-setting"  procedures,  such  as  Dljkstra's  algorithm[36] . 

Example  1.  Find  the  shortest  paths  from  node  1 to  all  other  nodes  in  the  undirected 
network  below. 


The  reader  is  encouraged  to  determine  the  shortest  paths  by  performing 
the  steps  indicated  in  Bellman's  and  Pape's  algorithms.  A crucial  computational 
consideration  in  both  Instances  is  the  length  of  list  T.  To  be  precise,  the 
number  of  elements  that  have  been  placed  on  the  list  will  determine  the  number 
of  executions  of  the  basic  computational  step,  the  most  costly  step  in  both 
procedures.  With  this  in  mind,  let  TB  be  the  list  of  nodes  for  the  Bellman 
algorithm  and  let  TP  be  the  list  of  nodes  for  the  Pape  algorithm.  To  avoid  con- 
fusion, we  performed  the  basic  step  in  ascending  order  of  k.  In  other  words,  we 
must  consider  arc  (1,  2)  before  arc  (1,  3),  and  so  on.  The  lists  are  given  below. 


TB  TP 


li 


I 


In  this  case,  Pape’s  algorithm  requires  almost  252  fewer  repetitions 
of  the  basic  computational  step  than  does  the  original  Hi*llman  algorithm.  In 
addition,  this  simple  example  is  indicative  of  more  general  problems  which 

Bellman's  algorithm  encounters  quite  frequently  and  which  Pape's  algorithm  is 
capable  of  avoiding.  Using  Bellman's  approach  In  example  1,  node  2 receives 
an  incorrect  minimal  distance  label  early  in  the  procedure  and  seven  other 
nodes  are  added  to  the  list  before  node  2 receives  Its  correct  minimal  label. 
Pape's  approach  adds  only  five  other  nodes  before  correcting  the  same  initial 
error.  In  general,  the  great  advantage  of  the  new  algorithm  is  that  errors  in 
minimal  distance  labels  are  corrected  as  soon  as  they  are  detected.  This  is 
accomplished  by  placing  the  node  with  the  corrected  label  at  the  top  rather  than 
bottom  of  the  list.  Pape  recommends  a "deque"  for  the  list  T and  he  discusses 
its  storage  as  well  as  computational  savings  (see  ^04] for  details).  A deque 
(or  double  ended  queue)  is  a linear  list  in  which  all  insertions  and  deletions 
are  made  at  the  ends  of  the  list. 

Problem  2 continues  to  generate  research  attention.  Golden  [ 66]  has 
studied  Problem  2 for  Euclidean  networks  only.  More  recently,  Denardo  and  Fox 
[33  ] have  introduced  a new  family  of  shortest  path  algorithms  based  on  buckets. 

A bucket  is  a list  of  nodes  whose  labels  fall  within  a given  range. 

Let  m (assumed  positive)  denote  the  length  of  the  shortest  arc  in  A. 

Then,  define  buckets  of  width  m such  Chat  bucket  p is  a list  of  nodes  1 whose 
temporary  labels  v(l)  fall  (currently)  in  the  interval. 

mp^v(l)<m(p  + l)»  p*l,  2,  ... 

In  Dijkstra's  algorithm  nodes  are  classified  either  as  permanently  or 
temporarily  labeled.  A permanently  labeled  node  is  one  with  a label  which  has 
been  shown  to  be  the  true  shortest  path  distance.  At  each  iteration,  Che 
algorithm  finds  the  node  with  Che  smallest  temporary  label  defined  by 
v(j)£  min  {v(l)  + d(i,j)  : i is  permanently  labeled  } and  makes  the  label 
permanent.  With  buckets,  we  can  replace  this  step  with  the  determination  of 
the  lowest-numbered  bucket  p*  Chat  contains  one  or  more  temporary  labels.  Suppose 
at  the  end  of  an  iteration  that  v(k)  is  Che  smallest  temporary  label.  Then,  by 
the  very  nature  of  phe  Dljkstra  algorithm,  all  nodes  1 such  that 
v(k)<^  v(i)  ^v(k)  + m 

must  be  permanently  labeled  (see  [ 33]  for  details).  Since  m is  the  length  of 
Che  shortest  arc,  it  is  Impossible  for  1 to  receive  a label  less  than  v(l) 
from  node  k or  any  teoiporarlly  labeled  node.  This  observation  can  result 
in  substantial  computational  savings. 
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We  note  that  deGhellinck  [32  ] has  had  encouraging  preliminary  computa- 
tional experience  imbedding  the  bucket  approach  to  shortest  paths  within  the 
out-of-kilter  algorithm  for  solving  transshipment  problems. 

The  research  that  we  have  been  discussing  is  primarily  at  the 
"implementation  level"  with  the  goal  of  developing  faster  and  faster  shortest 
path  algorithms.  Currently,  problems  with  thousands  of  nodes  are  being  solved  in 
fractions  of  a second.  Since  these  procedures  are  called  upon  routinely  by 
transportation  planners  in  so  many  applications,  this  development  is  worth 
following. 

The  other  four  shortest  path  problems  mentioned  earlier  have  also 
received  attention  in  recent  years.  Hart  et  al.t73],  Nemhauser  [ 98],  and  Golden 
and  Ball  fb7.]  discuss  the  application  of  a generalization  of  Dijkstra's  algorithm 
to  Problem  1.  Floyd's  algorithm  [ 50]  is  still  widely  cited  for  Problem  3.  As  an 
alternative,  we  can  repeat  Pape's  algorithm  from  each  node.  Dreyfus  [^0]  has 
proposed  an  algorithm  for  solving  Problem  4.  Kershenbaum  et  al.  [87]  also 
solve  this  problem  in  the  context  of  telephone  network  routing.  Regarding 
Problem  5,  Shier  [114  ],  ^15]  has  developed  an  extremely  effective  procedure  for 
finding  the  k shortest  paths  from  a given  node  to  all  other  nodes  in  a network. 
Dantzig  et  al.  [30]  have  recently  studied  decomposition  techniques  for  solving 
large  scale  shortest  path  problems. 

We  recommend  Dreyfus  [40 ],  Gilsinn  and  Wltzgall  [65],  and  Christofides 
[ 25]  as  general  sources  of  information  and  references  on  shortest  path  problems. 
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3.  Transshipment  Models 

In  making  shortest  path  trip  selections,  travelers  assume  that  their 
actions  do  not  effect  one  another.  Whenever  transportation  systems  are  operat- 
ing near  capacity,  however,  there  are  congestion  delays  and  this  assumption 
becomes  untenable.  In  these  situations,  realistic  models  should  account  for 
interactions  between  the  users.  Another  extension  to  the  shortest  path  model 
is  the  classical  transshipment  problem  in  which  goods  are  to  be  moved  simultane- 
ously from  several  sources  to  several  destinations,  for  example,  when  empty  rail 
freight  cars  are  to  be  transported  from  their  current  rail  yard  locations  to 
other  yards  where  the  cars  are  needed  to  provide  hauling  services  [12A  ] . In  this 
section,  we  study  some  of  these  important  extensions  to  the  shortest  path  model, 
concentrating  on  recent  algorithmic  contributions  for  solution  of  the  transship- 
ment problem. 


Modeling  Considerations 

To  set  notation,  we  cast  the  transshipment  problem  as  follows: 


minimize  I 
i 

subject  to  Z 

j 


E 

j 


ij 


'"ij  '‘iJ 


(i  - 1.  2, 


n) 


0 < Xj^j  i u^j  all  arcs  (i,j). 

In  this  formulation,  which  models  the  flow  of  a homogeneous  good,  be  it  a class 
of  passengers,  freight,  vehicles,  or  crew  personnel,  the  decision  variable  x^^ 
is  the  flow  of  the  good  on  arc  (1,J),  c^^  is  a given  per  unit  flew  cost,  and 
u^j  is  a given  upper  bound  (possibly  u^^  ■ + •)  for  flow  on  arc  (i,J).  The 
quantity  b^  is  the  known  supply  at  node  i,  a negative  value  being  interpreted 
as  a demand.  In  practice,  few  of  the  possible  arcs  in  this  network  model  will 
correspond  to  links  in  the  underlying  transportation  network,  with  three  to  five 
times  as  many  links  as  nodes  being  typical.  Accordingly,  we  assume  that  the  sum- 
mations and  indexing  in  the  model  are  restricted  to  links  (i,J)  and  (r,i)  of  the 
physical  network.  This  characteristic  of  network  sparsity  implies  that  a forward 
star  representation  of  the  data  is  an  attractive  storage  scheme. 

In  the  past  few  years,  remarkable  advances  have  been  made  on  these  problems. 
Using  contemporary  special-purpose  network  codes,  it  is  now  possible  to  obtain 
solutions  up  to  two  orders  of  magnitude  faster  than  by  solving  these  problems 
with  general-purpose  cosmerclal  linear  prograoming  packages.  Solving  the  trans- 
shipment problem  efficiently,  like  solving  shortest  path  problems  efficiently. 
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Is  Important  for  two  reasons.  First,  the  transshipment  problem  realistically 

models  a number  of  distribution  and  transportation  decision-making  situations. 

The  empty  freight  car  redistribution  problem  mentioned  previously  in  this  section 

and  the  aircraft  scheduling  model  introduced  in  an  earlier  section  are  but  two 

examples.  Moreover,  the  repeated  solution  of  this  model  is  often  embedded  within 

procedures  for  solving  more  complex  transportation  problems. 

As  an  illustration  of  this  second  point,  consider  a multi-fleet  routing 

version  of  the  aircraft  scheduling  model  with,  say,  707 's,  727's,  and  747's  as 

airplane  types  in  the  fleet.  The  constraints  of  the  transshipment  model  with 

superscripts  k « 1,  2,  or  3 differentiating  between  plane  types  on  all  variables 
k k k k 

x^j  and  on  all  data  b^,  c^j , and  u^j , models  the  routing  of  each  individual  plane 
type.  The  multi-fleet  model  includes  additional  "bundle"  constraints  of  the  form 


'ij 


+ 

*ij  *1J 


< 1 


on  all  service  links  (i,j),  for  example,  a potential  flight  leg  connecting  Boston 
at  8 AM  to  Atlanta  at  noon.  The  bundle  constraint  ensures  that  only  one  aircraft 
type,  if  any,  provides  this  service.  In  this  example,  we  would  require,  as  well, 
integral  values  for  the  decision  variables  x^j . 

The  more  general  version  of  this  multiconmodity  flow  problem  involves 
2 commodity  types,  each  subject  to  its  own  transshipment  constraints.  Bundle 
constraints  Imposed  upon  certain  arcs  of  the  network  with  bundle  capacities,  not 
necessarily  equal  to  one,  model  interactions  between  the  commodities.  The  model 
might  be  formulated  as  a linear  program  or  might  be  formulated  as  an  integer 
program  and  solved  via  branch  and  bound  using  linear  programming.  In  either 
case,  two  different  solution  strategies,  each  of  which  decomposes  the  problem 
into  a number  of  transshipment  models,  are  candidates  for  solving  the  linear 
program.  Price-directive  decomposition  places  a value  (or  price)  on  the 
bundle  capacity  of  each  arc  (i,J)  and  "charges"  for  use  of  this  capacity  in  a 
modified  (Lagrangian)  objective  function 


k 

I Z I (c^. 
k-1  i J 


to  be  minimized.  Resource-directive  decomposition  allocates  the  capacity  of 

each  bundle  arc  (i,j)  among  the  commodities,  i.e.,  imposes  additional  capacities 

It  ^ Ic  It  Ic 

y^j  on  the  flow  variables  such  that  0 i x^^j  5 min  (y^j  . Feasible  allocations 

k 

of  the  bundle  capacities  require  that  Z " K...  Both  algorithms  operate  by 

k-1  J 

V 

fixing  values  of  the  new  variables  X^j  or  y^j  and  discarding  the  bundle  constralnta 
so  that  the  problem  then  separates  into K independent  transshipment  models,  one 
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for  each  commodity  type.  Iteratively,  the  values  of  these  variables  are  read- 
justed until  an  optimal  solution  to  the  problem  Is  computed.  The  price-directive 
decomposition  approach,  Implemented  as  Dantzlg-Wolfe  decomposition,  has  been  very 
successful  In  several  recent  computational  studies.  Assad  [6  ],  [7  ],  Kennlngton 
[84],  [85],  Kennlngton  and  Shalaby  [86],  and  Swoveland  |120]  discuss  further 
details  of  these  algorithms,  describe  a number  of  applications,  and  report  on 
computational  experience. 

An  alternative  approach  to  modeling  of  multicommodity  flow  problems 
accounts  for  Interactions  between  the  connodlty  types  by  Incorporating  congestion 
effects  Into  the  objective  function.  The  transshipment  constraints  for  each  com- 
modity are  modeled  as  before,  the  bundle  constraint  Is  eliminated,  and  the  objec- 
tive function  is  replaced  by 

Minimize  f(x)  , 

where  f Is  a nonlinear  ftincdon  of  the  vector  x of  flow  variables  . The 
modeling  of  urban  traffic  flow  leads  to  an  important  class  of  problems  of  this 
type.  In  this  setting,  each  traveler,  or  group  of  travelers,  moving  between  an 
origin  node,  such  as  a suburban  housing  community,  and  a destination  node,  such 
as  a work  zone  In  the  central  business  district.  Is  Identified  as  a commodity. 

The  traffic  delays  on  any  street  (l,j)  of  the  network  might  depend  on  the  total 


flow  Z x^T,  on  that  street.  We  discuss  several  modeling  possibilities  for  this 
k-1 

problem  In  the  next  section.  Having  defined  the  delay  on  each  link,  we  might 
choose,  as  a normative  model,  to  minimize  total  delay  In  the  network  or  some 
function  of  delay  such  as  total  fuel  consumption.  These  scenarios  assume  that 
a central  authority  (e.g. , a city  planner)  directs  the  optimization  and  assigns 
routing  plans  to  all  travelers.  For  this  reason,  the  model  Is  often  referred  to 
as  system  optimization;  we  discuss  a decentralized  declslon-maklng  model  known 
as  user  optimization  In  Che  next  section.  As  we  shall  see,  user  optimization 
frequently  leads  to  Che  same  type  of  multicommodity  flow  model. 

We  postpone  discussing  solution  techniques  for  this  nonlinear  multicom- 
modity flow  problem  until  after  we  have  Introduced  Che  user  optimized  problem 
and  described  urban  traffic  flow  modeling  In  greater  detail  In  the  next  section. 
At  this  point,  we  merely  note  Chat  one  algorithmic  strategy  Is  to  linearize  the 
objective  function  and  repetitively  solve  transshipment  problems,  which  at 
times,  are  simple  shortest  path  problems. 
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Consequently,  solving  both  the  bundle  and  congestion  forms  of  the  multi- 
commodity flow  problem  requires  repetitive  and,  hence,  efficient  solution  of  the 
single  commodity  transshipment  model — a topic  that  we  consider  next. 

Solving  the  Transshipment  Problem  Efficiently 

Because  of  the  special  network  structure  of  the  transshipment  model, 
considerable  streamlining  is  possible  in  Implementing  the  simplex  method  to 
solve  the  problem.  Figure  3 shows  a typical  Iteration  of  the  simplex  method 
when  applied  to  a IS  node  transshipment  model. 


Adding  Arc  (8,  5)  to 
the  Basis 

Figure  3.  A Simplex  Iteration 

The  solid  arcs  in  this  figure,  these  arcs  correspond  to  the  variables 
x^j  in  the  linear  programnlng  basis.  Illustrate  a fundamental  and  well-known 
property  of  this  class  of  problems:  every  linear  programming  basis  corresponds 
to  a spanning  tree  in  the  underlying  network.  That  is,  when  arc  orientations 
are  ignored,  (1)  the  basic  arcs  contain  no  circuit,  and  (li)  any  nonbaslc  arc 
forms  a unique  circuit  with  the  basic  arcs.  The  darkened  arcs  in  Figure  3 are 
Che  basic  arcs  in  Che  circuit  formed  by  Che  nonbasic  arc  (8,  5). 
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In  Che  simplex  algorithm,  the  basis  is  updated  from  step  to  step  by 
introducing  a nonbasic  arc  to  replace  one  of  the  basic  arcs.  This  update  requires: 

(1)  determining  the  circuit  formed  by  the  basis  and  the  arc  being  intro- 
duced, which  we  call  the  pivot  circuit  (knowledge  of  the  circuit  and 
the  problem  data  determines  Che  arc  to  leave  Che  basis) , and 

(2)  recomputing  the  simplex  multipliers,  or  node  potentials,  that 
satisfy 

" ■ 'ij  - 'l  *3 

for  every  arc  (i,j)  in  the  new  basis. 

At  each  step,  any  nonbasic  arc  (i,j)  with  c^^  = c^j  ~ ^i  ^ ^ ^ becomes  a can- 

didate to  enter  the  basis.  The  method  of  selecting  from  these  candidate  arcs,  though 
something  of  an  art,  has  a profound  effect  upon  solution  time.  The  most  successful 
methods  choose  a subset  of  arcs,  varying  from  step  to  step,  and  introduce  Che  arc 
(p,q)  whose  reduced  cost  c^^  is  minimal  within  the  subset.  Mulvey  [97]  and  Bradley 
eC  al.  [20]  describe  several  mechanisms  for  implementing  this  strategy. 

Efficient  implementation  of  requirements  (1)  and  (2)  necessitates  the  care- 
ful storage  and  manipulation  of  data  describing  the  current  basis.  Since  these  issues 
have  been  so  successful  in  improving  algorithmic  performance,  and  since  they  illus- 
trate so  nicely  the  use  of  computer  science  techniques  in  transportation  appli- 
cations, we  describe  the  implementation  details  more  fully. 

In  Figure  3,  we  have  arbitrarily  "rooted"  the  basis  at  node  4.  Concep- 
tualizing a basis  in  this  manner  actually  facilitates  implementation.  In  order 
to  determine  the  pivot  circuit  formed  by  nonbasic  arc  (p,q),  we,  first,  may  find 
the  unique  paths  and  in  the  basis  connecting  these  nodes  to  the  root  of  the 
tree.  Those  arcs  lying  in  just  one  of  these  paths  together  with  (p,q)  form  the 
pivot  circuit.  Here,  it  is  convenient  to  record  the  predecessor  of  each  node  (the 
root  has  none) , which  is  the  first  node  encountered  when  traveling  from  this  node 
to  the  root,  i.e.,  the  next  higher  node  in  the  tree.  Then,  to  determine  the  paths 
Pp  and  P^  we  simply  trace  predecessors  to  the  root. 

Computing  the  paths  P^  and  P^  and  merging  them  to  find  the  pivot  circuit  is 

inefficient  in  most  circumstances,  but  particularly  when  the  circuit  is  small  and 

lies  deep  In  the  tree  so  that  the  paths  P and  P are  long  and  share  many  common 

P 9 

arcs.  An  alternative  is  to  move  upward  through  the  tree  from  nodes  p and  q,  one 
step  at  a time,  until  their  paths  meet.  There  are  several  ways  to  implement  this 
strategy.  For  example,  we  can  store  the  depth  of  each  node  in  the 
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tree;  the  depth  of  a node  is  the  number  of  arcs  In  the  path  Joining  the  node 
to  the  root.  In  Figure  3,  nodes  8 and  5 have  respective  depths  of  3 and  5. 
Because  the  paths  and  must  meet  at  the  same  depth,  we  may  find  the  pivot 
circuit  as  follows.  Move  from  the  deeper  of  the  nodes  p and  q,  say  p,  toward 
the  root  to  a node  that  is  as  deep  as  node  q.  Then,  move  towards  the  root  con- 
currently on  both  paths  P^  and  P^  until  the  point  where  the  paths  first  meet. 

An  alternative  is  to  store  the  number  of  successors  of  each  node  (i.e.,  the 
number  of  nodes  lying  below  it  in  the  tree).  We  then  move  from  the  node  with 
fewer  successors  (choose  arbitrarily  when  ties  arise)  towards  the  root  until 
we  again  encounter  the  same  node  on  both  paths  Pp  and  P^.  Both  methods  have 
been  implemented  successfully. 

Having  found  the  pivot  circuit  and  next  determined  the  outgoing  basic  arc 
by  one  of  these  techniques,  we  must  then  compute  the  simplex  multipliers  for 
the  new  basis.  Note  that  dropping  the  outgoing  arc,  arc  (12,  7)  in  Figure  3, 
splits  the  current  basic  tree  into  two  subtrees.  If  we  hold  the  simplex  multi- 
pliers for  all  nodes  in  one  of  these  subtrees  at  their  current  values,  then  to 

achieve  c..-tT.  + ti,  =0  for  all  arcs  in  the  new  basis  (and  c - tt  + tt  in 
ij  1 j pq  P q 

particular),  the  simplex  multipliers  for  every  node  in  the  other  subtree  must 
change  in  magnitude  by  Thus,  to  complete  this  step,  we  need  to  enumerate 

all  nodes  in  one  of  the  subtrees.  It  is,  of  course,  attractive  to  enumerate  the 
nodes  in  the  smaller  of  the  subtrees,  an  identification  that  is  simple  to  make 
when  the  number  of  successor  nodes  has  been  stored.  Sequencing  the  nodes 
properly,  and  maintaining  this  information  from  one  iteration  to  the  next, 
greatly  facilitates  enumerating  the  nodes  of  a subtree.  One  possibility  is  a 
sequence,  or  traversal,  that  "walks"  through  the  nodes  of  the  tree,  starting 
with  the  root,  from  top  to  bottom  and  left  to  right.  For  our  example,  this 
sequencing  would  read  4-9-15-8-1-12-10-7-11-5-3-14-2- 
6-13  before  and  4-9-  15  -8-5-  11  -7-3-  14  -1-12  -10  -2-6- 
13  after  the  basis  change.  These  traversals  satisfy  two  conditions:  (i)  the 
predecessor  of  each  node  appears  in  the  sequence  before  the  node  itself,  and 
(11)  directly  following  each  node.  In  sequence,  are  its  successors  in  the  tree 
(if  there  are  any). 

Now  suppose  that  we  wish  to  perform  the  basis  change  indicated  in  Figure  3. 
Deleting  the  outgoing  arc  (12,  7)  creates  two  subtrees.  To  enumerate  the  nodes 
below  node  7 in  one  subtree  and  update  the  simplex  multipliers,  we  extract  from 
the  traversal  the  sub-sequence  7-11-5-3-14.  Knowing  the  number  of 
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successors  of  node  7 deteimlnes  the  length  of  this  sub-sequence,  l.e.,  that 
it  terminates  at  node  lA.  This  information  can  also  be  obtained  from  the  pre- 
decessor data  since  2 is  the  first  node  after  node  7 in  the  sequence  whose 
predecessor  is  not  in  the  sub-sequence,  or  could  be  recorded  directly — for  each 
node,  store  the  last  node  in  the  sequence  that  is  one  of  its  descendents.  The 
choice  among  these  options  depends  upon  trade-offs  between  storage  and  computation 
time. 

We  have  now  seen  how  several  types  of  data  structures (predecessor , depth  or 
number  of  successors,  traversal)  might  reduce  computation  time  for  the  simplex 
method.  The  efficiency  of  the  algorithm  as  a whole,  however,  also  requires 
efficient  updating  of  these  structures  from  step  to  step.  The  change  in  the 
tree  is  rather  simple  conceptually.  "Holding"  the  tree,  with  the  Incoming  arc 
attached,  at  its  root,  we  "cut"  the  outgoing  arc.  The  tree  then  "falls"  into 
its  new  position  (see  Figure  3).  Note  that  the  subtree  below  the  cut  in  our 
example  appears  in  the  new  tree  with  the  path  from  node  3 to  the  cut  at  node  7 
reversed,  but  the  rest  of  the  subtree  remains  unchanged.  Exploiting  this 
observation  helps  in  updating  the  data  structures  efficiently. 

State-of-the-art  papers  by  Barr  et  al.  [lO]  and  by  Bradley  et  al.  1,20] 
describe  thoroughly  this  updating  process,  as  well  as  other  details  of  the 
algorithm  and  its  implementation.  These  papers  and  an  earlier  survey  by 
Magnanti  [94]  cite  and  review  a number  of  previous  contributions.  In  a related 
development,  Aashtlani  and  Magnanti  [ 2]  have  used  similar  data  structures  to 
reduce  the  computation  time  of  the  out-of-kilter  method  for  solving  the  trans- 
shipment problem. 

Theoretical  Bounds  on  Efficiency 

The  implementations  Just  described  have  proven  to  be  effective  in  solving 
numerous  problems;  they  lead  to  computation  times  of  about  8 seconds  on  a CDC  6600 
computer  for  solving  1500  node,  4300-5700  arc  problems.  Nevertheless,  Zadeh  (1261, 
(127 J has  constructed  arbitrarily  large  examples  requiring  a number  of  iterations 
that  is  exponential  in  N the  number  of  nodes.  These  examples  show  that  solution 
times  can  become  prohlbatlve  as  the  problem  parameter  N becomes  large.  Algorithms 
for  network  problems  like  the  transshipment  problem  are  said  to  be  "good"  if  the 
solution  time  for  any  example  is  bounded  by  a polynomial  in  N.  Several  researchers 
(Edmonds  and  Karp  [41],  Dlnlc  [37],  Karzanov  [83])  have  proposed  good  algorithms 
for  the  transshipment  problem,  or  special  cases.  One  result  of  their  effort  is  a 
novel  algorithm  for  the  maximal  flow  problem  whose  running  time  is  bounded  by  an 
order  N polynomial.  Even  [44]  reviews  the  algorithm  in  detail  and  Baratz[9] 
shows  that  this  run  time  bound  is  best  possible.  This  analysis  invites  further 
investigations  into  the  complexity  of  the  transshipment  model. 
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4.  Traffic  Equilibrium 

How  can  a subway  system,  an  expanded  major  artery,  a new  bridge,  one- 
way street  assignments,  priority  lanes,  and  other  policy  alternatives  available 
to  urban  planners  help  to  alleviate  congestion  in  our  cities?  More  specifically, 
how  would  users  respond  to  these  alternatives?  What  demands  would  they  impose 
upon  public  transportation  facilities,  what  routes  would  they  select  in  their 
travel  by  private  vehicles,  and  what  levels  of  congestion  would  the  system 
experience.  In  this  section,  we  consider  models  and  algorithms  for  predicting 
such  behavior. 

Ingredients 

The  flow  pattern  of  an  urban  transportation  network  depends,  to  a large 
extent,  on  relationships  between  demand  and  congestion: 

(1)  as  the  number  of  users  of  any  link  (arc)  of  the  transportation  net- 
work increases,  the  delay  time  (impedence)  along  that  arc  increases, 

and 

(2)  as  the  delay  times*  Increase,  the  demands  of  users  for  travel  decrease. 
The  models  that  we  consider  attempt  to  predict  flow  by  determining  when  the 
demand  "forces"  and  delay  time  "impedences"  equilibrate. 

In  practice,  the  most  successful  applications  of  urban  transportation 
modeling  have  been  limited  to  the  modeling  of  a single  transport  mode,  namely 
private  vehicles.  In  this  case,  the  delay  time  t^  along  an  arc  a is  frequently 
modeled  as 

tg(f)  = tpll  + a(f/Cg)^]  (1) 

Here  f is  the  total  flow  of  vehicles  on  the  arc,  c is  the  steady  state  capacity 
of  the  arc,  tQ  is  the  free  flow  time  and  a and  0 are  constants;  values  of 
a ■ 0.15  and  & » 4 are  typical  of  those  used  in  practice.  Branston  [2l]  describes 
a number  of  alternate  formulations  for  the  link  delay  curve  t (f)  and  makes 

a 

several  suggestions  concerning  the  proper  use  of  these  functions  in  practice. 

A weakness  of  delay  function  (1)  is  that  it  does  not  account  for  the 
fact  that  the  delay  along  a link  is  often  a function  of  flows  on  other  links  in 
the  network.  For  example,  the  delay  along  a link  feeding  into  a busy  intersection 
might  depend  upon  flow  on  other  links  feeding  that  intersection.  Furthermore, 
since  two-way  streets  are  modeled  as  two  (directed)  arcs  with  opposite  orienta- 
tions, the  delay  in  one  direction  is  often  a function  of  the  traffic  in  the  other 
direction  due,  in  part,  to  left  hand  turns. 


* Any  user  cost  may  be  used  in  place  of  delay  time  throughout  this  discussion. 
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The  demand  component  of  urban  transportation  modeling  usually  concen- 
trates on  origin  and  destination  points,  0-D  pairs,  for  travel.  Households, 
businesses,  and  other  end  points  for  travel  are  aggregated  into  zones  that  are 
represented  by  nodes  in  the  transportation  network.  Other  nodes  in  the  network, 
such  as  intersections,  are  transshipment  points  for  vehicle  travel. 

Most  systems  currently  available  for  predicting  urban  traffic  flow,  such 
as  the  UMTA  (Urban  Mass  Transit  Authority)  Transportation  Planning  System,  generate 
demands  by  considering  trip  production  and  attraction  factors  such  as  income  and 
parking  availability  in  the  zones,  and  travel  time  between  the  zones.  Having  fixed 
interzonal  demands  by  this  trip  distribution  phase,  a trip  assignment  procedure, 
such  as  the  equilibrium  model  that  we  discuss  in  the  next  subsection,  predicts  the 
route  choice  that  users  make  in  order  to  meet  their  travel  demands. 


Because  the  traffic  patterns  generated  in  this  second  phase  may  provide 
new  estimates  of  travel  times  between  the  zones,  adjustments  to  demands  may  be 
called  for.  The  ultimate  prediction  of  urban  flow  would  be  obtained  by  iterat- 
ing between  the  trip  distribution  and  trip  assignment  phases  in  some  way,  either 
formally  or  heuristically . 

This  iteration  can  be  automated  by  using  demand  functions  in  the  equilib- 
rium model  that  depend  upon  travel  times  with  respect  to  prevailing  network  con- 
gestion. We  might  model  demand  between  each  0-D  pair  1 as  a function  D^(u^), 
possibly  linear,  of  the  shortest  travel  time  u^  between  that  pair.  More 
generally,  demand  could  be  expressed  as  D^(u),  a function  of  the  vector  u of 

shortest  travel  times  Uj  between  all  0-D  pairs  j • 1,  2 n.  This  extended 

formulation  permits  broader  modeling  capabilities,  such  as  incorporating  destina- 
tion choice.  Suppose,  for  example,  that  0-D  pairs  1 and  2 represent  travel  from 
a given  zone  to  each  of  two  shopping  districts.  Dial's  [35]  extended  "logit 
model"  with 


Dj^(u)  - d 


0Ui 
r^e  1 


0Ui  . 0U5 

J-  + r e Z 


+ r^e 


, DjCu)  - d 


0U2 


0Ul  . 0U9 

tj^e  + r2e  ^ 


(2) 


where  d is  the  total  number  of  shopping  trips  to  be  made  and  r^^  and  r2  are 
attraction  factors  for  the  two  shopping  districts,  permits  the  traffic  assignment 
procedure  make  destination  choices  between  the  shipping  centers. 

When  time  dependent  demand  models  are  used,  the  second  phase  procedure 
simultaneously  determines  traffic  distribution  and  traffic  assignment  with 
factors  such  as  zonal  income  and  parking  availability  held  fixed.  These  models 
are  short  range  planning  tools.  Longer  range  analysis,  for  Instance,  studies  on 
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the  impact  of  new  transportation  facilities  on  urban  development,  would  require 
data  relating  variations  in  the  "fixed"  factors  to  demand. 

For  further  discussion  of  transportation  demand  models,  the  reader  might 
consult  Domenenich  and  McFadden's  monograph  [39].  Nguyen  [lO]] , [102  ] considers 
a problem  related  to  our  discussion,  namely  estimating  0-D  zonal  trips  from 
observed  link  flows. 


An  Equilibrium  Model 

In  his  seminal  paper  [122],  Wardrop  posited  two  principles  for  determining 
the  distribution  of  traffic  in  an  urban  network.  The  first  of  these,  which  is 
the  basis  for  most  urban  transportation  planning  models,  is  a fundamental  behav- 
ioral assumption  about  user  objectives  in  traveling  between  a given  0-D  pair: 

"the  journey  time  on  all  routes  actually  used  are  equal,  and  less 
than  those  that  would  be  experienced  by  a single  vehicle  on  any 
unused  route." 

The  following  mathematical  model  of  equilibrium  captures  the  principle: 


Tp(h)  > u. 


all  p e P, 


T (h)  = if  h > 0,  p E V (i“l,2 n) 


Z h - D (u) 
peP^  P 


(3) 

W 

(5) 


where  T (h)  * Z {t  (h) : arc  a belongs  to  path  p}. 

P acA  ^ 

This  model  is  usually  referred  to  as  a user  equilibrium  and  contrasts  with  the 
system  equilibrium  model  introduced  in  the  last  section  which  corresponds  to 
Wardrop 's  second  principle. 

In  the  formulation  (3)-(5),  A denotes  the  arcs  of  the  network,  i“l,2, 
...,n  denotes  the  0-D  pairs,  and  P^  denotes  a set  of  paths  joining  0-D  pair  1. 
Usually  consists  of  all  paths  joining  0-D  pair  i,  although  it  may  contain  a 
subset  of  these  paths  (e.g.,  only  those  that  the  user  perceives).  The  term  hp 
is  the  number  of  vehicles  using  path  p;  u^^,  D^(u)  and  t^(h)  are  the  shortest 
travel  times,  demands  and  link  delays  introduced  previously  Tp(h),  the  sum  of 
travel  times  along  all  links  a belonging  to  path  p,  gives  the  total  travel  time 
on  path  p. 

Condition  (3)  states  that  the  travel  time  on  any  path  joining  an  0-D  pair 
must  be  at  least  as  large  as  the  shortest  travel  time.  Condition  (4)  states  that 
a user  travels  (denoted  by  hp  > 0)  only  on  paths  giving  the  shortest  travel  time 
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ir^ 


between  any  0-D  pair.  Equality  (5)  implies  that  all  demand  for  travel  between 
any  0-D  pair  is  satisfied  by  flow  along  paths  Joining  that  0-D  pair. 

As  Rosenthal  pOT]  has  pointed  out,  this  formulation  is,  in  a sense,  a 
continuous  approximation  to  Wardrop's  principle,  since  in  practice  vehicles  are 
indivisible  and  the  path  flow  variables  should  be  integral.  Welntraub  [12j  has 
studied  relationships  between  the  continuous  and  Integral  formulations.  Since 
the  integer  model  has  not  been  Implemented,  we  shall  confine  our  discussion  to 
the  continuous  model. 

Several  features  of  this  model  are  worth  noting.  First,  we  observe  that 
the  link  delay  function  for  each  arc  a is  written  in  terms  of  the  entire  flow 
pattern  in  the  network  and  not  merely  in  terms  of  the  total  flow 


f^  a 2{hp:  arc  a belongs  to  path  p} 


(6) 


on  that  arc.  Thus  the  model  can,  in  principle,  incorporate  the  link  inter- 
actions mentioned  in  the  previous  section.  Second,  because  each  demand  function 
D^(u)  is  expressed  in  terms  of  the  shortest  path  distances  between  all  0-D  pairs, 
the  model  has  the  potential  to  provide  for  destination  choice  modeling  and  other 
extensions. 

Finally,  we  should  recognize  that  the  formulation  (3)- (5)  encompasses 
multimodal  distribution.  To  Illustrate  this  point,  let  us  consider  modeling  of 
buses  and  autos  as  two  alternate  means  of  transportation.  We  envisage  two  net- 
works, one  for  each  mode  (possibly  copies  of  the  same  road  network)  identifying 
a mode  with  Its  network.  Each  path  set  then  represents  0-D  transport  by  one 
mode.  The  demand  and  link  delay  function  will  embody  interactions  between  the 
modes.  In  this  simple  setting,  0-D  Indices  1 and  2 might  correspond  to  bus  and 
auto  travel  between  the  same  physical  0-D  pair;  the  logit  model  (2)  is  one  pos- 


sibility for  expressing  the  demands  D^(u}  and  D2(u}  for  the  two  modes.  The 
interpretation  of  (2)  is  much  the  same  as  before;  d is  the  total  demand  between 
the  0-D  pair  and  r^  and  r2  are  attraction  factors  for  the  modes. 

Recently,  Florian  [43}  and  Abdulaal  and  LeBlanc  [ 3}  have  proposed  two 
mode  equilibrium  models  along  these  lines  and  suggested  algorithms  to  compute 
a solution.  Note  that  this  type  of  model  is  capable  of  computing  traffic  dis- 
tribution, modal  split,  and  traffic  assignment  simultaneously,  in  contrast  to 
the  one  pass  sequential  approach  of  the  widely-used  UMTA  Transportation  Planning 
system.  For  related  work  see  Bruynooghe  [22],  Florian,  Nguyen,  and  Ferland 
[49],  Evans  [43]  and,  particularly,  Florian  and  Nguyen  [ 48] • 
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Although,  as  we  have  seen,  the  equilibrium  formulation  permits 

richness  In  modeling,  calibrating  the  link  delay  functions  and  demand  functions 
and  computing  equilibria  for  the  most  general  formulation  remains  an  unattalned 
objective.  One  comforting  feature  of  this  modeling  approach  Is  that  very  mild 
restrictions  on  the  problem  data  (t  (0)  > 0,  t (h*)  t t (h)  whenever  h'  i h, 
t^(h)  continuous,  and  D^(u)  continuous  and  bounded  from  above)  guarantee  that 
an  equlllbrum  exists.  Aashtlanl  [ 1 ] recently  established  this  fact  using 
results  from  nonlinear  complementarity  theory. 

The  use  of  the  equilibrium  model  as  a planning  tool  has,  to  date,  been 
limited  to  single  mode  private  vehicle  applications  In  which  (1)  the  volume 
delay  on  each  link  depends  only  on  the  total  flow  f^  on  that  link  as,  for 
example.  In  (1),  (11)  the  demand  between  each  0-D  pair  depends  solely  upon  Che 
shortest  travel  time  between  that  origin  and  destination,  l.e.,  demand  Is  given 
by  Dj^(u^),  and  (ill)  D^(u^)  is  a decreasing  function.  The  key  to  analyzing 
this  situation  Is  an  observation  made  by  Beckman,  McGuire,  and  Winston  [ ll], 
that  the  Kuhn-Tucker  conditions  to  the  following  optimization  problem  (in 


variables  h and  d.) 

P i 


Minimize  1 
acA 


(t)  dx 


n ^*1 
- I / £l 


(y)  dy 


subject  to  £ 
peP, 


1,  2,  . . .,  n) 


h i 0.  d^  i 0 peP^,  (1-1,  2, 


are  equivalent  to  the  equilibrium  conditions  (3)- (5)  when  the  Kuhn-Tucker  multi- 
plier A.  for  Che  1^^  equality  constraint  is  Identified  with  the  shortest  travel 
time  Uj^  between  0-D  pair  1,  if  d^^  > 0.  In  this  optimization  model  g^(y)  ■ (y) 

Is  the  Inverse  of  the  demand  function  (the  model  Includes  the  special,  but 
Important)  case  of  constant  demands  D^(u^)  by  setting  g^(y)  * 0). 

The  importance  of  modeling  the  equilibrium  problem  as  an  equivalent 
minimization  problem  Is  that  the  objective  function  of  (7)  is  convex  whenever 
t^(T)  and  gj^(y)  fulfill  the  practical  assumptions  of  being,  respectively,  non- 
increasing  and  nondecreasing.  Consequently,  methods  from  convex  programming  can 
be  applied  computationally. 

Computing  An  Equilibria 

Of  Che  several  proposals  that  have  been  made  for  computing  equilibria  by 
solving  (7),  see  Nguyen  L99],  the  most  widely  used  is  the  Frank-Wolfe  algorithm. 
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This  method  solves  nonlinear  programs  with  linear  constraints,  l.e., 
min  {f(x):  Ax  “ b,  x i O),  by  repeated  linearization  of  the  objective  function. 
Given  any  feasible  solution  x-*  to  the  problem,  the  method  finds  a solution  y 
to  the  linearized  problem 

min  {Vf  (x^)  • y:  Ay  ■ b,  y i 0}  (8) 

where  Vf(x^)  is  the  gradient  of  f evaluated  at  . It  then  solves  the  one- 
dlmenslonal  search  problem  of  minimizing  f In  the  line  segment  joining  x^  and 
y,  obtaining  a new  solution  x^^^.  The  method  Iterates  over  j » 0,  1,  2,  . . . 
starting  with  an  arbitrary  Initial  feasible  solution  x^. 

This  algorithm  is  particularly  well  suited  for  solving  problem  (7). 
Consider,  first,  the  fixed  demand  model  In  which  d^  Is  a constant  and  g^(y)  ■ 0 


for  1 =■  1,  2,  . 


n 

I Z 

n.  Any  linear  objective  function  p C h is 


minimized  subject  to  the  constraints  of  (7)  by  setting  l^q  “ *^1 

satisfies  Cq^  ■ min  {Cp  : pcP^};  that  is,  the  demand  is  met  by  any  minimum  cost 

path.  Some  algebraic  manipulations  reveal  that  the  coefficient  Cp  of  hp 

obtained  by  linearizing  (7)  about  any  vector  (h  of  given  path  flows  is 

■j  P 

simply  the  sum  of  t (f  '^)  along  arcs  a belonging  to  path  p.  Consequently,  the 

A d 

linearized  problem  (8)  reduces  to  a sequence  of  shortest  path  problems,  one  for 
each  0~D  pair,  with  the  prevailing  link  delays  as  arc  costs.  These  shortest 
path  problems  can  be  solved  efficiently  by  the  techniques  described  In  an 
earlier  section  of  this  paper. 

Two  points  about  the  algorithm  are  worth  noting.  First,  there  Is  no 
need  to  enumerate  all  paths  In  each  path  set  P.  prior  to  the  analysis.  The 
algorithm  generates  them  as  needed.  Second,  only  the  total  flow  t on  each 

O 

arc  needs  to  be  maintained  from  step  to  step.  Once  the  one-dlmenslonal  line 
search  has  been  performed  at  each  step,  the  shortest  path  solutions  can  be  dis- 
carded. Exploiting  this  fact  leads  to  substantial  reductions  In  storage 
requirements. 

The  variable  demand  version  of  (7)  Is  solved  In  much  the  same  way. 

When  linearized,  the  objective  function  coefficient  of  for  problem  (7) 
becomes  u^  ■ g^fd^^).  The  solution  to  the  subproblem  (8)  then  depends  upon 
both  the  values  of  Cq^as  defined  above,  and  the  current  shortest  path  distances 
u^.  A solution  Is: 
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where  is  any  known  upper  bound  on  the  demand  between  0-D  pair  i.  Nguyen  [lOO] 
describes  further  details  about  this  algorithm  and  discusses  other  methods  for 
solving  for  an  equilibrium  with  variable  demands. 

Several  researchers  have  contributed  ideas  related  to  this  algorithm. 

The  excellent  survey  [A7 ] contains  additional  references  and  provides  his- 
torical perspective  concerning  this  and  other  algorithms  for  solving  the  minimiza- 
tion problem  (7).  We  should  emphasize  that  problem  (7)  is  just  one  manifestation 
of  the  congestion  formulation  of  the  multicommodity  flow  problems  discussed  in  the 
last  section.  Any  algorithm  for  solving  (7)  usually  applies  to  this  broad  generic 
set  of  models. 

An  alternative  to  casting  the  equilibrium  problem  in  equivalent  convex 
minimization  form  (7)  is  to  view  the  model  as  a nonlinear  complementarity  problem. 
The  model  then  can  be  studied  from  the  viewpoint  of  this  theory  (Aashtlani  [ 1 ], 
Hall  [72  ])  or  the  viewpoint  of  fixed  point  theory  (Kuhn  [91  ], Kuhn  and  Cullum  [92]). 

This  approach  has  the  advantage  of  applying  to  the  general  equilibrium  formulation, 
but  the  disadvantage,  to  date,  of  requiring  much  greater  computer  time  and  storage 
than  the  minimization  approach. 

Computational  Experience 

To  test  the  validity  of  equilibrium  modeling  as  a predictive  tool, 

Florlan  and  Nguyen  [46]  applied  model  (7)  to  data  from  the  city  of  Winnipeg. 

They  assumed  fixed  travel  demands,  as  generated  from  a previous  study,  and  used 
an  alternate  to  (1)  for  modeling  link  delays.  The  model  predicted  flow  on  high 
volume  links  quite  well,  but  did  not  perform  as  well  on  links  with  observed 
volumes  in  the  range  of  0-300  vehicles  per  hour.  Their  findings  show,  as  might 
be  expected,  that  the  predictions  of  route  travel  times  were  better  than  those 
of  link  travel  times.  They  concluded  that  "the  results  are  encouraging  and 
demonstrate  the  suitability  of  the  method  for  planning  purposes." 

In  this  study  the  Frank-Wolfc  algorithm,  equipped  with  a Dijkstra-type 
shortest  path  routine,  required  lS-18  iterations  and  about  700  CPU  seconds  on 
a CDC  Cyber  74.  The  convex  simplex  method  solved  the  same  problem  in  about 
500  CPU  seconds,  but  required  more  storage.  The  network  contained  1319  links. 
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In  another  study.  Hern  [ 78]  considered  a 9386  link,  3027  node  network 
of  Washington,  D.C.  The  Frank-Wolfe  algorithm,  as  Implemented  In  the  TRAFFIC 
computer  code  , required  221  CPU  seconds  per  iteration  on  an  IBM  360/91. 


I 


28 


5.  Vehicle  Routing 

Like  many  operational  Issues  In  transportation  planning,  the  vehicle  routing 
problem  Is  encountered  routinely  and  repeatedly  In  business  and  industry.  The 
basic  problem  is  one  of  designing  a set  of  vehicle  routes  of  minimal  total 
distance  leaving  from,  and  eventually  returning  to,  a central  depot,  which 
satisfies  capacity  constraints  and  meets  customer  demands.  Demands  occur  at 
points  or  nodes  in  the  transportation  network  and  may  be  deterministic  or 
probabilistic  In  nature.  Generally,  there  are  enormous  amounts  of  detailed 
data  and  a far  larger  number  of  feasible  sets  of  routes  to  consider.  As  a 
result,  only  small  problems  can  be  solved  for  optimal  solutions;  otherwise,  we 
must  reconcile  ourselves  to  the  fact  that  heuristic  solutions  (hopefully  near- 
optimal)  must  suffice.  In  addition,  there  are  a host  of  Inter-related  aspects 
of  the  vehicle  routing  problem  Including  the  number  and  location  of  depots  and 
demand  points,  the  capacity  of  vehicles  and  makeup  of  fleet,  frequency  of 
service,  and  other  geographical  considerations.  The  computational  complexity 
of  this  problem  and  the  fact  that  this  type  of  problem  is  often  solved  every 
day  underscores  the  need  for  powerful  and  efficient  solution  techniques. 

Deterministic  Setting 

First,  we  focus  on  the  case  wliere  demands  are  deterministic.  Examples 
of  this  problem  include  municipal  waste  collection  [13],  fuel  oil  delivery  [55], 
newspaper  distribution  [69],  and  routing  of  school  buses  [14],  Notice  that  in 
some  examples  pick-ups  are  made;  in  others  deliveries  are  made.  As  long  as  only 
one  of  the  two  operations  is  performed  throughout,  the  distinction  is  not  important. 

There  are  many  heuristic  techniques  which  have  been  proposed  for  this  class 
of  problems.  We  concentrate,  in  this  section,  on  one  approach  which  has  been  suc- 
cessful in  solving  large  problems  (more  than  100  nodes).  Gillette  and  Miller  [64] 
and  Orloff  [103]  discuss  alternative  heuristic  strategies. 

The  algorithm  we  describe  is  an  efficient  implementation  of  the  Clarke- 
Wrlght  savings  method  [27].  Suppose  we  let  node  1 denote  the  central  depot  and 
d(l,j)  be  the  distance  from  node  1 to  node  J.  If  every  two  demand  points  i and 
j are  supplied  individually  by  two  vehicles  from  the  central  depot,  then  total 
distance  traveled  is  2 d(l,l)  + 2 d(l,J).  However,  if  both  points  are  served 
by  a single  vehicle  then  the  combined  route  results  in  a savings  in  travel 
distance  of 

{2  d(l,i)  + 2 d(l,j)}  - {d(l,l)  + d(l,j)} 

- d(l,i)  + d(l,j)  - d(l,j). 
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A similar  savings  occurs  whenever  Che  endpoints  of  two  vehicle  routes  are 
Joined  to  form  a single  route.  In  Che  Clarke-Wrlght  algorithm,  we  proceed  as 
follows: 

Step  1.  Evaluate  all  potential  savings 

S(i,j)  - d(l,i)  + d(l,j)  - d(i,J)  for  i,j  i 1. 

Step  2.  Order  all  feasible  savings  from  largest  to  smallest. 

Step  3.  Select  the  node  pair  (1,J)  with  the  greatest  positive  feasible  savings. 
Link  nodes  i and  j on  a single  tour. 

Step  4.  Eliminate  infeasible  savings  and  return  to  step  2,  until  there  are  no 
remaining  positive  feasible  savings. 

When  we  say  a savings  S(i,j)  is  feasible  we  mean  that  linking  nodes  i and  j does 
not  cause  Che  violation  of  any  constraint  (tour  integrity,  vehicle  capacity, 
maximum  route  time,  and  so  forth).  After  each  linking  of  nodes,  a number  of 
savings  become  infeasible  (see  [69]  for  details).  For  example,  an  intermediate 
tour  l-2'l-3-l  implies  that,  for  all  k,  savings  S(l,k)  are  infeasible  since 
otherwise  we  would  not  preserve  Che  tour. 

This  algorithm  can  be  coded  efficiently  by  taking  advantage  of  two 
important  observations: 

1.  Step  1 can  be  very  costly  both  computationally  and  in  terms  of  storage 
requirements.  For  instance,  a 600-node  problem  requires  360,000  storage 
locations  for  Inputs  (distances  d(l,J)  above  Che  diagonal  and  savings 
S(1,J)  below  the  diagonal  for  an  undirected  network  with  symmetric  dis- 
tances). 

2.  At  each  step  of  the  algorithm,  we  must  determine  the  maximum  feasible 
savings. 

These  observations  can  be  exploited  by  using  special  data  structures 
and  list  processing  techniques.  First,  rather  than  consider  an  entire  matrix 
of  pairwise  linkings,  we  can  focus  on  the  most  promising  linkings  only. 

Instead  of  storing  the  network  topology  in  a matrix,  we  would  record  for  each 
potential  arc  its  origin  node,  its  destination  node,  and  its  length,  l.e. , use 
a ladder  representation.  From  this  Information  we  then  calculate  the  savings. 

We  restrict  the  entries  In  this  list  to  reduce  the  number  of  savings  con- 
sidered to  under  10,000  for  a 600-node  problem.  One  procedure  for  accomplish- 
ing this  reduction  Is  to  consider  linking  node  1 to  any  node  j within  a 
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distance  of  r units  of  i.  Golden,  Magnanti,  and  Nguyen  [69]  document  an 
alternative  approach  involving  a rectangular  grid.  An  extremely  convenient 
and  efficient  method  for  finding  the  best  feasible  savings  is  to  partially 
order  the  savings  in  a heap  structure  and  update  the  structure  from  step  to 
step  (see  [69])«  These  ideas  lead  to  an  implementation  of  the  Clarke-Wrlght 
algorithm  which  is  between  one  and  two  orders  of  magnitude  faster  than  the 
traditional  implementation. 

Stochastic  Setting 

Now,  we  consider  the  more  complex  problem  of  constructing  a fixed  set 
of  routes  when  demands  are  probabilistic.  Such  a problem  would  arise  when 
dally  deliveries  of  fuel  oil  are  being  made  to  automotive  service  stations 
and,  although  each  route  is  fixed  in  advance,  the  demand  on  a particular  day 
is  stochastic.  For  simplicity,  let  us  assume  that  the  demand  at  each  node  i, 
denoted  by  d^,  can  be  modeled  by  a Poisson  distribution  with  mean  The 

discussion  here  follows  Stewart  [119]  and  Golden  and  Stewart  [70]. 

We  say  that  a primary  error  has  occurred  if  a vehicle  cannot  satisfy 
the  demands  of  the  customers  on  the  route  to  which  it  has  been  assigned.  This 
situation  has  various  penalty  costs  associated  with  it.  Clearly,  one  objec- 
tive is  to  minimize  the  probability  of  a primary  error.  The  stochastic 
vehicle  routing  problem  can  then  be  formulated  as  determining  a fixed  set  of 
routes  to: 

Minimize  (1)  expected  total  travel  distance 
subject  to  (2)  meeting  customer  demands; 

(3)  not  exceeding  vehicle  capacity; 

(4)  Prob  {primary  error  on  a route}  1 a. 

We  can  solve  the  problem  heurlstlcally  in  such  a way  that  we  take 
advantage  of  the  efficient  Clarke-Wrlght  implementation  discussed  in  connec- 
tion with  deterministic  demands.  Suppose  a route  contains  nodes  n^^,  n2, 
n^  and  has  total  demand  and  total  expected  demand 

x*d  +d  +...+ 

"l  "2 

E(x)  ■ X + X + ...  + X , respectively. 

"l  "2  "hi 

Then,  by  appealing  to  the  Central  Limit  Theorem  (see  [70]  for  details)  we  can 
approximate  the  Poisson  distribution  for  total  demand  by  a Normal  distribution 
with 
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W-X  +A  +...+A  and 
“l  "2  “k 

o ■ /in 

If  we  assume  chat  all  vehicles  have  the  same  functional  capacity  c,  then  the 
probability  of  a primary  error  Is  given  by 

Prob  (x  > c)  - Prob  {z  I ~ 

/in 

using  the  Nomuvl  approximation  where  z Is  a unit  normal  variate. 

Our  solution  strategy  will  be  to  replace  the  functional  capacity  of  a 
vehicle  with  a reduced  "artificial"  capacity  and  to  replace  the  stochastic 
demand  at  each  node  with  a deterministic  "artificial"  demand  In  such  a way 
Chat  Che  deterministic  model  may  be  used. 

Let  y denote  the  artificial  capacity  of  a vehicle  and  A^  the  artificial 
demand  at  node  1.  If  vehicles  are  loaded  to  their  artificial  capacities,  Chen, 
solving 

Prob  iz  - — } = o , 

/r 

we  obtain,  after  some  algebra,  Che  value  of  y which  Insures  Chat  constraint 
(4)  Is  satisfied.  That  Is, 

- ^i-g  ~ ^l-g 

^ ~ 2 

where  is  defined  by  Prob  {z  - *2.-0^  ■ 1 - a.  For  Instance,  If  c • 100 

and  a * .10,  Chen  > 1.28  and  y ~ 87.9.  Using  an  Integral  artificial 
capacity  of  87  units,  this  gives  a safety  stock  of  13  units  as  a cushion 
against  the  occurrence  of  primary  errors.  Fixed  routes  are  constructed  on 
the  basis  of  the  artificial  capacity  and  artificial  demands,  as  In  the  deter- 
ministic case.  Sensitivity  analysis  Is  recommended  for  differing  values  of  a. 
Golden  and  Stewart  [70]  have  Implemented  this  solution  strategy  and  presented 
computational  results.  We  point  out  Chat  the  procedure  applies  to  many  prob- 
ability distributions  other  than  the  Poisson. 

A potential  transportation  application  Involves  the  design  of  an  effec- 
tive subscription  bus  service  for  large  employment  centers.  A reliable  dally 
bus  service  would  make  pickups  at  pre-speclfied  bus  stops  each  morning,  travel 
to  the  employment  center,  and  make  drops  at  the  bus  stops  In  the  evening.  There 
might  be  a fixed  monthly  fee  for  service.  Since  a subscriber  need  not  show  on 
a particular  day  because  of  Illness,  vacation,  or  special  plans,  we  have  a 
probabilistic  vehicle  routing  problem. 
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6.  Network  Design 

There  are  a number  of  options  in  the  design  of  a transportation  system, 
ranging  from  operational  alternatives  such  as  one-way  street  assignments  to  the 
physical  improvement  of  existing  facilities  or  construction  of  new  facilities. 
We  shall  view  any  of  these  options  as  a network  synthesis  and  refer  to  any 
alternative  in  terms  of  a network  construction. 

A general  taxonomy  of  design  problems  divides  Into  (1)  arc  construction 
In  which  certain  arcs  (e.g.,  roadways  or  rallbeds)  are  added  to  the  network,  or 
not,  and  (11)  facility  location  models  in  which  nodes  representing  warehouses, 
depots  and  the  like  are  "opened,"  or  not.  In  each  case,  the  objective  Is  to 
determine  economic  tradeoffs  between  the  cost  of  construction  and  the  savings 
In  routing  cost  that  it  provides. 

After  briefly  reviewing  several  models  for  assessing  these  tradeoffs, 
we  consider.  In  this  section,  recent  algorithmic  advances  for  this  class  of 
problems.  Our  discussion  focuses  on  three  different  approaches — branch  and 
bound.  Benders  decomposition,  and  heuristic  procedures. 


Design  Models 

For  k a 1,  2,  ...,  K let  denote  the  flow  requirement  between  nodes 
and  tj^  of  a given  network;  let  c^.  denote  the  per  unit  routing  cost  on  arc 


(i,j)  for  "commodity  k"  goods,  and  let  f^^j 
Ing  arc  (l,j).  Then,  a rather  general  network  design  model  Is: 


denote  the  fixed  cost  of  construct- 


Minimlze  ZEE  c^  xj^  + Z E f y.  . 

k 1 j J 1 j J 


subject  to 


E x^.  - E x^. 
j ij  , rl 


f 'k  “ ^ ■ *k 

< --^k  “ ‘ ■ ^k 
I 0 otherwise 


ij  ^Ij 


[ j ®ij  ® 

(x,y)  e S 


X > 0 

*ij  - ^ 


^ij 


all  1 and  J 
• 0 or  1 all  i and  J . 


(9) 

(10) 
(11) 
(12) 
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0)  and  cannot  exceed  the  capacity  of 


In  this  formulation,  is  the  flow  on  arc  (i,j)  of  goods  being  transported 
from  node  Sj^  to  node  tj^.  The  binary  variables  y^j  indicate  whether  or  not  an 
arc  is  constructed.  Indexing  conventions  for  this  model  are  similar  to  those 
we  have  used  for  the  transshipment  model  in  section  3. 

Equations  (9)  are  the  usual  transshipment  constraints.  The  "bundle" 
Inequalities  (10)  specify  that  the  total  flow  on  arc  (i,j)  must  be  zero  if  that 
arc  is  not  constructed  (i.e.,  y^ 

the  arc  if  it  is  constructed.  Inequality  (11)  is  a budget  constraint  stating 
that  total  construction  cost  is  limited  by  a budget  B.  In  this  expression, 

Cjj  is  a cost  function,  which  need  not  equal  ^ t related  to  the  building  of  arc 
(i,j).  The  set  S encompasses  further  side  constraints,  possibly  expressed  as  linear 
inequalities  in  the  vectors  x = and/or  y ■ These  might  include  pre- 

cedence relations  (e.g.,  construct  arc  (i,j)  only  if  arc  (i',j')  is  constructed), 
multiple  choice  relations  (e.g.,  choose  at  most  (at  least,  exactly)  two  of  some 
subset  of  arcs),  limitations  on  resources  shared  by  several  arcs,  or  prior 
specifications  of  arcs  already  constructed  such  as  y^^  * 1.  When  prior  speci- 
fications have  been  made,  the  model  is  often  referred  to  as  a network  improve- 
ment model. 

Most  papers  written  about  network  design  consider  specializations  of 
this  model.  We  shall  use  the  following  terminology  for  these  problem  variants: 

Uncapacitated  Design — every  K. . ^ ^ r,  so  that  constraints  (10)  impose 

ij  k 

no  (capacity)  restrictions  on  flow  when  y^^  “ 1. 

Fixed  Charge  Design — the  budget  constraint  is  eliminated: 
tradeoffs  between  fixed  charge  and  routing  costs  are 
investigated . 

Budget  Design — the  fixed  charge  term  E eliminated. 

Optimal  routing  is  sought  within  a fixed  construction  budget. 

We  should  note  that  although  our  formulation  Includes  only  one  construction  level 
for  each  arc,  multiple  capacity  levels  can  be  modeled  by  parallel  links.  Also, 
as  an  alternative  to  the  bundle  constraints  (11),  we  might  incorporate  nonlinear 
congestion  costs  in  the  objective  function.  This  type  of  formulation  is  attractive 
for  continuous  models  where  Incremental  improvements  are  being  made  to  an  existing 
network.  Dantzlg  et  al.  [29],  [30]  and  Harvey  and  Robinson  [74]  have  considered 
such  a model  for  budget  design.  They  apply  a Lagrange  multiplier  to  the  budget 
constraint,  use  the  Frank-Wolfe  algorithm  (see  section  3)  for  any  fixed  multiplier 
value,  and  Iterate  to  find  the  optimal  value  for  the  multiplier. 
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A facility  location  model  similar  to  this  formulation  of  the  network, 
design  problem  Is: 

Minimize  L E c,.  x. . + E f y 
1 j IJ  1 ^ ^ 

subject  to  E = dj  (13) 


E 

j 


X 


ij 


(lU) 


E y^  < B (15) 


(x,y)  e S 

1 0 all  i and  j 

= 0 or  1. 

In  this  case,  the  underlying  network  is  bipartite.  There  are  m possible  loca- 
tions for  (production)  facilities.  The  variable  y^  equals  1 if  site  i is 
selected  and  is  zero  otherwise.  denotes  the  capacity  of  a site  if  it  is 

selected,  d.  denotes  the  demand  at  destination  node  J and  f.  is  a fixed  cost 
J ^ 

for  selecting  site  i.  As  before,  (15)  represents  a budget  constraint  and  the 
set  S incorporates  various  side  conditions.  When  the  fixed  charge  term  of  the 
objective  function  is  omitted,  each  e^  = 1,  and  B = P is  a limit  on  the  number 
of  facilities  that  can  be  selected,  the  model  is  called  a P-median  problem. 

As  demonstrated  by  Wong  [i25  ] , the  location  model  can  be  cast  as  a net- 
work design  problem  by  adding  an  artificial  node  and  an  arc  joining  this  node 


to  every  source  node  1.  Let  r^  denote  the  flow  requirement  from  this  new  node 
to  destination  J and  associate  f^,  K^,  and  e^  with  the  arc  joining  the  new  node 
with  source  node  i. 
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Johnson  et  al.  [81]  have  shovm  that  the  network  design  budget  problem, 
even  with  each  e^j  *1,  is  NP-complete.  That  is,  it  can  be  solved  by  an 
algorithm  that  is  polynomial  in  problem  input  (number  of  nodes,  size  of 
budget)  if  and  only  if  any  of  a number  of  other  notoriously  difficult  problems 
(including  the  traveling  salesman  problem,  and  the  multicommodity  flow  problem 
in  integers)  can  be  solved  similarly.  This  result  provides  some  theoretical 
insight  concerning  the  difficulty  of  design  problems. 

Branch  and  Bound 


Several  researchers  have  proposed  branch  and  bound  algorithms  for  solv- 
ing network  design  and  facility  location  models  (Boyce  et  al.  [ 1^,  Chrlstofldes 
and  Brooker  [26],  Dionne  and  Florian  [38],  Leblanc  [93]  Scott  [ill]  Steenbrink 
[118],  Efroymson  and  Ray  [42],  Davis  and  Ray  [ 3l]»  Hoang  [79]  and  Geoffrion 
and  McBride  [63],  among  others).  To  illustrate  the  nature  of  th.is  work,  we 
describe  two  of  the  more  recent  contributions  from  this  list. 

Several  years  ago,  Hoang  [ 79]  suggested  an  enumerative  algorithm  for 
budget  design  problems  with  unit  flow  requirements  between  every  pair  of  nodes. 


In  this  Instance,  the  optimal  routing,  given  any  network  configuration,  is  via 
shortest  distance  paths  Joining  each  origin-destination  pair.  Consequently, 
the  objective  function  cost,  denoted  F(y),  is  determined  completely  by  the 
network  configuration,  i.e.,  by  the  choice  of  0-1  values  for  the  components 


y . of  the  vector  y.  At  any  node  P in  the  branch  and  bound  enumeration  tree, 

F <*  o . 

certain  arcs  A are  fixed  as  constructed  (either  y^j  = 1 or  y^^  ■*  0). 

As  a bounding  mechanism  for  his  algorithm,  Hoang  noted  that  any  solution 

y with  the  arcs  in  A fixed  at  these  same  values  satisfies  the  inequality 


F(y)  i F(y^)  + E y 1..  (y^) 


(16) 


- P 

where  y^j  ■ 1 - y^j.  In  this  expression,  y denotes  a solution  with  y^j  * 1 
for  every  arc  not  in  A^  and  = y^j  for  every  arc  (l,j)  e A^.  Thus  F(y'^)  is 
the  best  possible  shortest  route  solution  with  the  given  fixed  values  of  arcs 
in  A^.  I^j  (y^)  denotes  the  increment  to  the  shortest  route  cost  from  node  1 
to  node  J when  arc  (l,j)  is  deleted  from  the  network  defined  by  y^. 

Expression  (16)  has  the  following  interpretation.  If  arc  (1,J)  is 
deleted  (set  y^j  ■ 0 and  y^j  ■ 1)  from  the  network  defined  by  y^,  then  the 
cost  for  shipping  the  unit  of  demand  between  these  nodes  must  Increase  by  at 
least  I^j(y^)  in  the  solution  y.  It  might  increase  by  more  because  other  arcs 
are  being  deleted  as  well.  Hence,  the  right-hand  side  of  (16)  is  a lower  bound 
on  the  routing  cost  F(y)  of  solution  y. 


3b 


To  compute  lower  bounds  quickly,  Hoang  suggests  relaxing  the  integer 

requirement  on  and  minimizing  the  right-hand  side  of  (16)  subject  to  the 

budgetary  constraint  L L e^^  y^^  1 a and  0 ^ y. ^ ^ 1,  y^j  =•  y^j  for  arcs 

(i,j)  c A^.  A solution^y*  to  this  continuous  knapsack,  problem  has  at  most 

one  fractional  component  and,  in  light  of  (16),  gives  a lower  bound 

F(y)  i F(y^)  + 1 y*..  1..  (y^)  that  applies  to  every  node  below  node 

fi  i ) ^ A^  ^3  ^3 

P in  the  branch^adn^ bound  enumeration  tree.  Using  this  lower  bound  as  a 
fathoming  mechanism  and  branching  from  node  P on  (i.e.,  next  fixing)  the 
fractionally  valued  variable  in  the  continuous  knapsack  solution,  Hoang  imple- 
ments his  algorithm  in  the  framework  of  a straightforward  branch  and  bound 
algorithm. 

Dionne  and  Florian  [ 38]  have  noted  several  ways  to  improve  this  algo- 

p 

rithm.  First,  in  computing  I.,  (y  ) it  is  not  necessary  to  resolve  from  scratch 
for  shortest  paths  between  all  nodes.  Specialized  algorithms  are  available  for 
recomputing  shortest  path  distances  when  one  arc  has  been  deleted  from  a net- 
work. Second,  they  suggest  branching  on  the  variable  y^j « (i.j)  ^ with  high- 
est incremental  improvement  per  unit  of  budget,  i.e.,  that  arc  (i,j)  maximizing 
(y^)/e^^.  These  modifications  lead  to  marked  improvements  in  the  algorithm. 

A typical  example  with  a 65%  budget  level,  i.e.,  B = 0.65  Z Z with  20  nodes, 

i j ^3 

and  with  30  arcs,  requires  20  seconds  to  solve  on  a CDC  Cyber  74  computer  with 
their  algorithm,  while  Hoang's  algorithm,  after  500  seconds,  is  not  able  to 
determine  that  the  current  best  solution  is  optimal.  The  authors  note,  however, 
that  this  branch  and  bound  approach  is  probably  limited  to  medium-sized  networks 
like  this,  and  that  computation  time  seems  to  grow  exponentially  with  a decrease 
in  budget  level.  For  example,  the  algorithm  requires  288  seconds  to  solve  the 
same  problem  with  a 50%  budget  level. 

In  another  development  Geoff rion  and  McBride  1 63]  consider  Lagrangean 
relaxation  embedded  within  a branch  and  bound  approach  to  the  fixed  charge 
facility  location  problem.  To  adhere  to  their  notation,  let  us  assume  that 
dj  > 0 for  all  j and  let  us  substitute  z^^  s formulation.  They 

assume  that  the  side  conditions  of  S are  expressed  as  a system  of  linear 
inequalities  Ay  + Bz  £ b and  attach  a vector  of  Lagrange  multipliers  u to  these 
constraints  and  a vector  of  Lagrange  multipliers  X to  the  constraint  (13)  to 
form  a Lagrangean  relaxation: 
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L(X,u) 


Minimize  Z Z c'  . . z 


y.z 

subject  to: 


i i 


+ Zf^y^  + ZX  (Z^Z^^-l)  + u(b  - Ay  - Bz) 

i j ^ 


ij  ij 


Z dj  z,.  . K,  y, 
1 1 


0 ^ z 


ij 


all  i 

all  i and  j 
0 or  1 all  i 


(17) 

(18) 
(19) 


Here, 


ij  -ij  "z 

There  are  several  reasons  for  considering  this  type  of  relaxation.  The 

well-known  "weak  duality"  property  of  Lagrangean  duality  demonstrates  that 

L(X,w)  is  less  than  or  equal  to  the  optimal  value  v of  the  location  problem 

for  any  value  of  X and  any  u - 0.  Thus,  the  Lagrangean  relaxation  provides  a 

lower  bound  chat  can  be  utilized  as  a bounding  mechanism  in  branch  and  bound. 

This  particular  Lagrangean  relcixation  is  attractive  because  it  is  so  easy  to 

solve.  It  separates  into  independent  subproblems,  one  for  each  i.  For  each  i, 

either  y^  ••  0 and  z^^j  =•  0 for  all  j or,  y^  » 1 and  the  computation  of  optimal 

z^j  reduces  to  a continuous  knapsack  problem.  Moreover,  the  value  L(X,’,i)  of  Che 

Lagrangean  relaxation  when  X and  u are  set  equal  to  optimal  dual  values  of  the 

linear  programming  relaxation  provides  a sharper  (larger)  lower  bound  on  v 

than  does  the  optimal  value  of  Che  linear  programming  relaxation. 

Computing  the  tightest  Lagrangean  lower  bound,  i.e.,  finding 

L H max{L(X,y)  : u k 0}  yields  insight  into  the  structure  of  location  models. 

A y 

Due  to  tne  equivalence  between  dualization  and  convexif ication  (see  Geoffrion 
t 57]  and  Magnanti  et  al.  L^^])»  L also  equals  the  optimal  objective  value  of  Che 
problem: 

Minimize  Z Z c..  z..  + Zf^y. 

z,y  i j ^ ^ i ^ 

subject  to  Z z ■ 1 all  j 

4 


A y+  B z f b 


and  (x»y)  e Convex  Hull  of  solutions  to  (17),  (18)  and  (19). 

As  Geoffrion  and  McBride  point  out,  Che  convex  hull  constraints  of  this  formula- 
tion are  equivalent  to  the  linear  Inequalities  (17),  (18),  and  0 i y^^  - 1, 
together  with  the  inequalities  z^j  i y^  for  all  i and  j.  This  result  gives  some 
theoretical  justification  for  appending  Che  constraints  z^^  i y^  for  all  1 and  j 
Co  Che  location  model.  Moreover,  Geoffrion  and  McBride  conjecture  chat  "linear 
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programming  technology  will  advance  to  the  point"  that  solving  for  L in  this 
linear  inequality  form  may  be  "preferable  to  involving  Lagrangean  relaxation." 

Computational  experience  with  problems  ranging  from  7 possible  facili- 
ties and  122  transportation  links  to  25  possible  facilities  and  1250  links 
indicates  that  Lagrangean  relaxation  provides  much  tighter  lower  bounds  on  v 
than  does  linear  programming  relaxation.  A branch  and  bound  algorithm  equipped 
with  Lagrangean  relaxation  solved  these  problems  in  from  3.4  to  112.9  seconds 
on  an  IBM  370/158  computer  compared  with  a range  of  from  6.8  to  greater  than 
300  seconds  for  a branch  and  bound  algorithm  equipped  with  conventional 
Driebeek-Tomlin  penalties. 


Benders  Decomposition 


When  applied  to  ae'work  design  problems,  Benders  decomposition  proceeds 
iteratively  by  choosing  a tentative  network  configuration  (i.e.,  setting  values 
for  the  integer  variables  7^j)»  solving  for  the  optimal  routing  on  this  network, 
and  using  the  solution  to  the  routing  problem  in  order  to  redefine  the  network 
configuration.  Figure  4 illustrates  this  last  step  for  an  uncapacitated  fixed 
cost  design  problem  in  which  one  unit  of  a good  is  to  be  sent  from  node  1 to 
node  6.  In  this  instance,  the  routing  problem  reduces  to  a shortest  path  compu- 
tation between  nodes  1 and  6. 


Tti  « 0 100  1^6  ■ 100 


Figure  4.  A Step  of  Benders  Decomposition 

Tlie  dark  arcs  in  the  figure  are  members  of  the  current  network  conflgura-  ' 

i 

tlon;  the  dashed  arcs  are  candidates  for  inclusion  in  the  optimal  design.  With  ' 

respect  to  the  routing  costs  showr  next  to  each  arc,  the  node  numbers  n are  ^ 

optimal  dual  variables  for  the  linear  programming  routing  problem.  The  dual 
variables  indicate  that  introducing  arc  (2,  3}  into  the  current  design  reduces 
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the  routing  cost  from  node  1 to  node  3 from  ^3  ■ 50  to  ^<2  ^23  “20  + 10  for 

a savings  of  50  - 30  - 20  units.  Similarly,  introducing  arc  (4,  5)  reduces 
the  routing  cost  from  node  4 to  node  6 by  15  - (1T4  + C45)  » 90  - (70  + 10)  ■ 

10  units.  Since  either  of  these  arcs  might,  but  need  not,  become  part  of  the 
shortest  route  path  from  node  1 to  node  6 in  the  optimal  network  design,  the 
optimal  routing  cost  R is  constrained  by 

R > 100  - 20  yzz  - 10  yu5  • ^20) 

That  is,  at  best,  the  current  routing  cost  which  is  100  units  will  be  reduced 
by  the  full  savings  of  introducing  arc  (2,  3)  (i.e.,  setting  yzz  “ D the 
full  savings  of  introducing  arc  (4,  5)  (i.e.,  setting  yi,5  * 1). 

Constraints  like  (20),  which  are  known  as  Benders  cuts,  are  by-products 
of  the  optimal  routing  calculation  for  any  tentative  network  configuration. 
Benders  algorithm  computes  the  new  configuration  at  each  step  by  minimizing 
the  fixed  charge  design  cost  F plus  the  routing  cost  bound  R subject  to  the 
Benders  cuts  (20)  generated  by  every  previous  tentative  configuration.  This 
minimization,  called  the  Benders  master  problem,  is  a mixed  integer  program  in 
the  integer  variables  y^^  and  single  continuous  variable  R.  At  each  step,  one 
new  constraint,  a Benders  cut,  is  added  to  the  master  proLlam.  Note  that  since 
R becomes  a lower  bound  on  the  routing  cost,  the  minimum  cost  v of  any  design 
is  bounded  by  the  optimal  value  F*  + R*  of  the  master  problem,  i.e.,  v 1 F*+  R*. 
Also,  every  solution  y ■ y to  the  master  problem  determines  a network  and  the 
combined  fixed  and  optimal  routing  cost  on  this  network  is  an  upper  bound  on 
V.  These  two  bounds  permit  early  termination  of  the  algorithm  with  an  assess- 
ment of  degree  of  suboptimality . 

When  applied  to  more  complicated  design  problems  such  as  facility  location 
problems,  and  even  to  general  mixed  integer  programs.  Benders  decomposition 
operates  in  much  the  same  way  and  has  similar  interpretations. 

In  their  highly  regarded  paper  [61  ] (see  also  [62])  concerning  a facility 
location  model  with  shipments  from  plants  to  customers  through  intermediate 
distribution  centers,  Geoffrion  and  Graves  report  on  the  most  successful  docu- 
mented implementation  of  Benders  decomposition  to  date.  They  solve  problems 
with  from  249  to  513  binary  variables.  From  0 to  30  of  these  correspond  to 
distribution  centers  to  be  opened,  or  not;  the  remainder  are  0-1  variables 
indicating  whether  or  not  a distribution  center  serves  a customer.  The  compu- 
tations required  no  more  than  seven  Benders  iterations  to  reach  within  0.20% 
of  the  minimum  cost  design  value,  and  required  from  16.7  to  191  seconds  of 
execution  time  on  an  IBM  360/91. 
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Their  paper  is  rich  in  its  description  of  implementation  both  in  terms  of 
management/model  interaction  and  computer  programming  of  Benders  algorithm. 

The  authors  note,  for  example,  that  solving  the  Benders  master  problem  to  com- 
pletion at  each  iteration  may  be  an  unnecessary  computational  burden.  Rather, 
they  search  for  a solution  at  each  step  that  merely  increases  the  lower  bound 
by  at  least  a given  constant  c > 0.  They  also  note  that  modeling  constraints 

like  Z X i K y . , where  K ■ Z d , in  the  uncapacitated  facility  location 
j ij  1 i i j J 

model  (see  (14))  as  x^^  - y^  for  all  i and  j leads  to  much  better  algorithmic 
performance  despite  the  fact  that  the  latter  representation  requires  many  more 
constraints.  We  comment  further  on  this  observation  in  the  next  subsection. 

Even  for  a given  model  representation,  it  is  possible  to  accelerate 
Benders  decomposition  by  generating  "strong  cuts"  at  each  iteration.  Referring 
to  Figure  4 will  help  illustrate  this  point.  Note  that  the  shortest  path 
I distance  from  node  3 to  node  6 using  all  arcs  that  are  candidates  for  Inclusion 

in  the  optimal  network  design  is  60  units.  Since  the  distance  to  node  2 in  the 
current  design,  as  specified  by  the  dark  arcs,  is  2C  units  and  the  current 
shortest  path  cost  is  100  units,  introducing  arc  (2,  3)  whose  routing  cost  is 
10  units  can  save  no  more  than  100  - (60  + 20  + 10)  ■ 10  units,  and  not  the  20 
units  computed  earlier.  Consequently,  a valid  Benders  cut  is 

R 1 100  - 10  y23  - 10  yu5  . (21) 

Note  that  this  cut  is  stronger  in  the  sense  that  it  provides  a tighter  lower 
bound  on  R,  than  (20) ; the  right-hand  side  of  (21)  is  as  large  as  that  of  (20) 
for  all  0-1  values  of  the  decision  variables  y^^ , and  exceeds  the  right-hand 
' side  of  (20)  whenever  y23  ■ 1. 

The  opportunity  to  generate  strong  cuts,  like  (21),  is  made  possible 
^ because  of  degeneracy  in  the  shortest  path  linear  program,  or  equivalently 

, because  of  multiple  optimal  solutions  to  its  dual.  In  this  example,  ■ 0, 

7r2  ■ 20,  ff3  - 40,  - 70,  tts  ■ 90,  and  irg  ■ 100  is  an  alternate  optimal  dual 

solution  to  that  shown  in  Figure  4.  Computing  a Benders  cut  as  before,  but 
using  these  dual  values  leads  to  the  stronger  cut  (21).  Because  network  problems 
are  renowned  for  their  degeneracy,  considerations  of  this  nature  are  attractive 
for  a number  of  transportation  applications. 

Magnantl  and  Wong  [96  ] describe  this  strong  cut  methodology  in  greater 

detail.  They  show  how  to  generate  pareto-optimal  cuts,  i.e.,  cuts  with  the 

i 

[ property  that  no  other  is  stronger,  for  arbitrary  mixed  integer  programs  by 
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solving  a linear  program  to  choose  fiom  among  optimal  dual  solutions.  Special- 
izations of  this  approach  lead  to  a pareto-optimal  cut  methodology  for  P-median 
problems  that  avoids  explicit  linear  programming  computations.  Computational 
experience  on  a variety  of  P-median  problems  (up  to  33  nodes)  shows  that  Benders 
algorithm  equipped  with  strong  cuts  finds  solutions  known  to  be  within  102  of 
optimality  in  ten  or  fewer  iterations.  The  standard  implementation  usually  pro- 
vides no  better  solutions  within  25  iterations  and  solutions  102  farther  from 
optimality  within  ten  iterations.  The  authors  obtain  similar  experience  with 
strong  cuts  for  uncapacitated  fixed  charge  design  problems,  though  in  this  case 
the  error  bounds  are  generally  not  as  tight. 


Heuristics 

Several  researchers  (Billheimer  and  Cray  [16],  Scott  [ill],  Stairs  Lll7]» 
Steenbrink  [11^,  .and  Dionne  and  Florian  [38]  among  others)  have  proposed 
heuristic  procedures  for  solving  network  design  problems.  These  are  generally 
of  three  types — add,  delete,  or  interchange.  The  add  heuristics  start  with 
some  feasible  design  and  add  arcs,  one  at  a time,  choosing  at  each  stage  the 
arc  chat  gives  Che  greatest  decrease  in  cost,  or  some  surrogate  measure  of 
cost.  The  delete  heuristics  are  similar,  but  start  with  an  initial  design 
containing  all  candidate  arcs,  and  delete  arcs  one  at  a time.  Starting  with 
some  initial  design,  the  interchange  heuristics  add  and/or  delete  an  arc  at 
each  step  until  no  further  improvement  in  cost  is  possible. 

Recently,  Dionne  and  Florian  [38]  have  reported  impressive  computational 
experience  with  a new  type  of  heuristic  for  budget  design  problems  with,  for 
convenience,  one  unit  of  demand  between  every  two  nodes  in  the  network.  They 
use  their  branch  and  bound  algorithm  described  earlier  in  this  section  with  the 

p 

following  modification.  In  place  of  the  term  Ijj(y  ) in  the  lower  bound 

P 

expression  (16),  they  use  a term  I(y  } which  is  the  increment  to  shortest  route 
costs  between  every  pair  of  nodes,  not  Just  1 and  J,  when  arc  (i,J)  is  deleted 
from  Che  network.  This  algorithm  is  a heuristic  because,  as  examples  show,  the 
new  lower  bound  need  not  be  valid. 

The  authors  have  tested  this  algorithm  on  problems  ranging  from  7 nodes 
and  16  candidate  arcs  to  29  nodes  and  54  candidate  arcs , and  they  have  compared 
its  performance  with  an  add  heuristic,  with  Scott's  [lllj  combined  delete  and 
exchange  heuristic,  and  with  a variant  of  this  algorithm.  With  but  one  excep- 
tion, In  which  Che  relative  error  between  the  heuristic  and  optimal  solutions 
was  0.032,  Che  new  heuristic  found  Che  optimal  solution  to  every  problem.  The 


r 


42 


relative  errors  ranged  from  about  1%  to  7%  for  the  add  heuristic  and  were  less 
than  1%  in  all  cases  for  the  other  heuristics.  The  new  heuristic  required  from 
0.05  to  8.47  seconds  of  computation  time  on  a GDC  Cyber  74  computer. 

Cornuejols  et  al.  [28],  who  cite  a number  of  references  on  heuristics  for 
facility  location  models,  have  initiated  a new  line  of  investigation.  They 
consider  uncapacitated  fixed  charge  facility  location  models  with  the  P-median 

1 


constraint  1 y^  i P 


Let  us  consider  the  special  case  where  f^  ■ 0,  d^ 


and  1 0 for  all  i and  j.  These  assumptions  are  not  essential,  but  are 
merely  convenient  for  our  exposition.  To  conform  with  this  paper,  we  assume 
that  the  problem  has  been  formulated  in  maximization  form. 


The  authors  derive  the  following  results.  If  v is  Che  optimal  value  to 
Che  problem  and  v^  is  the  value  of  Che  solution  determined  by  the  add  heuristic 
Chen 


I 


e 


(22) 


where  e is  Che  base  of  the  natural  logarithm.  The  clever  proof  of  this  result 
uses  the  fact  that  the  value  v^  of  Che  Lagrangean  dual  problem  formed  by  dual- 
izing with  respect  to  the  demand  constraints  I x^j  “ 1 is  equivalent  to  the 

i 

optimal  value  of  the  linear  programming  version  of  the  problem  (i.e.,  0 i y^  5 1) 
when  Che  constraints  x^^  i y^  for  all  i and  j replace  (14).  The  analysis  of  Che 
Lagrangean  dual  shows  that  (22)  is  valid  when  v^  replaces  v.  Since  v®  £ v S vD, 
we  also  have 


Zf  --1  < (P  -l)P  < i 

vD  p « 


(23) 


which  indicates  by  how  much  Che  value  of  Che  linear  programming  relaxation  of 
Che  problem  can  deviate  from  the  value  of  Che  problem  itself. 

Not  only  do  Che  authors  show,  by  examples,  that  these  bounds  are  best  pos- 
sible, they  also  construct  examples  to  show  that  ocher  heuristics  (delete, 
exchange,  dynamic  programming)  cannot  provide  relative  error  bounds  as  tight  as 
(22).  Moreover,  they  show  chat  Che  relative  error  of  Che  optimal  value  of  Che 
"weak"  linear  progranmlng  relaxation  based  upon  Che  constraints  (14) , rather 
chan  Che  "strong"  formulation  based  upon  Che  constraints  x^^  i y^  for  all  1 
and  j , need  not  be  as  tight  as  (23)  (the  relative  error  can  be  mada  as  close  to 
1 as  desired  by  judicious  choice  of  data).  This  fact  helps  to  explain  modeling 
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experience  mentioned  previously,  namely  that  Che  strong  linear  programming 
formulation  is  preferred  to  the  weak  formulation. 

In  computational  experiments  with  problems  containing  as  many  as  164 
potential  locations  for  facilities,  Che  authors  used  Che  add  heuristic  followed 
by  a subgradienc  algorithm  applied  to  the  Lagrangean  dual  to  obtain  an  Improved 
feasible  solution  and  upper  bound.  This  algorithm  determined  Che  optimal  solu- 
tion to  almost  every  problem  chat  was  tested  and  required  no  more  chan  30  seconds 
of  computation  time  on  an  IBM  370/168  computer. 
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7 . Open  Research  Problems 

In  this  section,  we  enumerate  several  areas  where  future  work  may 

prove  fruitful: 

* Aggregation  of  Network  Data.  In  most  transportation  studies,  underlying 
data  Is  aggregated  In  order  to  obtain  data  sets  that  are  manageable  from 
both  the  viewpoints  of  data  collection  and  limitations  on  algorithmic 
capabilities.  What  are  the  "besi."  procedures  for  aggregating  data,  how 
are  aggregate  solutions  to  be  disaggregated  for  the  actual  problem  setting, 
and  what  degree  of  suboptimality  does  the  aggregation/disaggregation 
process  entail?  chan  et  al.  [23],  Geoffrion  [58],  [59],  Kuhn  [91],  and 
Kuhn  and  Cullum  [ 92  ] have  initiated  research  to  answer  these  questions. 

* Shortest  Paths  from  an  Origin  to  a Few  Destinations.  There  are  very  effi- 
cient algorithms  for  finding  shortest  paths  from  an  origin  to  one  other 
node  and  to  all  other  nodes.  In  practice,  however,  it  might  be  desired 
that  several  nodes  serve  as  destinations.  Can  this  problem  be  approached 
directly?  In  particular,  when  distances  are  Euclidean  can  the  geometry 

be  exploited  to  obtain  more  efficient  algorithms? 

* Multicommodity  Flow.  Can  list  processing  structures,  like  those  recently 
developed  for  transshipment  problems  (see  section  3),  be  extended  effi- 
ciently to  solve  multicommodity  flow  problems  with  bundle  constraints? 

Are  there  effective  heuristic  techniques  for  large  multicommodity  flow 
problems? 

* Normative  Models  of  Urban  Traffic  Flow.  Suppose  that  a system  optimized 

traffic  equilibrium  (e.g.,  fuel  minimization)  is  desired,  but  that  users 
act  in  accordance  with  Wardrop's  first  principle  of  individual  cost,  or 
time,  minimization.  What  available  policy  alternatives,  e.g.,  one-way 
street  assignments,  prohibited  turns,  or  road  tolls,  guarantee  that  a user 
equilibrium  would  coincide  with  the  desired  system  equilibrium?  Rosenthal 
[108]  has  shown  that  road  tolls,  alone,  will  not  suffice.  To  perform  analysis 
of  this  nature  might  require  solution  methodologies  which  Incorporate  pre- 
dictive models  as  part  of  the  underlying  constraint  structure  of  normative  I 

optimization.  Theoretical  tools  for  performing  sensitivity  analysis  on  j 

equilibrium,  solutions  (see  Hall  [7l])  would  undoubtedly  aid  this  effort,  and 

would  be  of  independent  Interest  to  transportation  planners.  | 
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* Predicting  Origin  and  Destination  Choice  of  Urban  Traffic.  When  the  traf- 
fic equilibrium  model  Includes  destination  choice  possibilities  (see  section 
4)  there  Is  no  known  equivalent  convex  optimization  problem  for  predicting 
traffic  flow.  Origin  choice  models  that  would  be  capable  of  predicting 
homeowner  location  decisions  for  long  range  urban  planning  encounter  the 
same  difficulty.  The  ability  to  compute  solutions  to  either  variant  of 
this  model,  or  to  a combined  origin  and  destination  choice  model,  would 
greatly  extend  urban  traffic  planning  capabilities. 

* The  Stochastic  Vehicle  Routing  Problem.  There  has  been  some  recent  work 
on  this  problem  as  described  In  this  paper,  but  much  more  research  needs 
to  be  done  and  real  applications  need  to  be  litvestlgated. 

* The  Vehicle  Routing  Problem  with  Window  Constraints.  Timing  restrictions 
become  a component  of  the  vehicle  routing  model  In  the  event  that  some 
customers  impose  delivery  deadlines  and  earliest  delivery  time  constraints, 
thereby  Imposing  a "window"  or  "interval"  of  time  during  which  a delivery 
or  pickup  must  be  made.  Biles  and  Bradford  [ 15]  and  Russell  [lO^  have 
recently  looked  at  these  sequencing  restrictions,  but  more  efficient 
algorithms  should  be  within  reach. 

* Tight  Lower  Bounds  for  Vehicle  Routing  Problems.  Can  a relaxation  approach 
similar  In  nature  to  the  work  of  Held  and  Karp  [76  1,  [77]  be  of  value  In 
deriving  tight  lower  bounds  on  the  minimum  distance  solution  to  the  vehicle 
routing  problem.  This  would  allow  one  to  evaluate  the  effectiveness  of 
various  heuristic  approaches. 

* Airline  Crew  Scheduling.  An  airline  has  a set  of  m flights  each  of  which 
requires  one  of  a set  of  N crews.  Crews  must  be  allocated  to  flights  in 
order  to  minimize  operating  costs  (see  Arabeyre  et  el.  [5]  for  details). 
Powerful  heuristics  need  to  be  developed  and  analyzed  for  this  Important 
class  of  problems. 

* Worst-Case  Analysis  of  Heuristic  Algorithms.  For  a given  network  optimiza- 
tion problem  %diere  the  only  known  efficient  solution  procedures  are  heuristics, 
how  badly  might  the  heuristic  solution  deviate  from  the  optimal  solution? 
Rosenkrantz  et  al.  [106]  and  Cornuejols  et  al.  [28]  provide  examples  of  this 
kind  of  analysis. 

* Probabilistic  Analysis  of  Network  Algorithms.  See  Karp  [g?]  for  a number  of 
open  problems  of  a more  theoretical  nature  relating  to  probabilistic,  as 
opposed  to  worst-case,  analysis  of  heuristics. 
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* The  Loading  Problem.  Let  x^j  and  y^j  denote  the  flow  of  goods  and  vehicles 

on  the  arcs  (i,j)  of  the  same  network.  Given  a linear  or  nonlinear  objec- 
tive function  of  x ■ y “ ^^ij^  minimized,  add  to  the 

separate  transshipment  constraints  for  goods  and  vehicles  loading  constraints 
of  the  form  x^^  - K for  each  arc  (i,j).  The  interpretation  is  that  goods 
are  to  be  loaded  on  vehicles  each  with  capacity  K.  Applications,  among 
others,  are  the  loading  of  passengers  on  planes,  the  loading  of  cargo  on 
trains,  and  the  assignment  of  railcars  to  rail  engines.  Straightforward 
modeling  extensions  encompassing  multicommodity  goods  and  multiple  vehicle 
types  might  be  expected  in  practice.  What  is  an  efficient  solution  tech- 
nique for  this  class  of  problems? 

* Freight  Flow  Management.  Modeling  of  freight  involves  a number  of  Issues. 

In  rail  applications,  decisions  must  be  made  concerning  the  assignment  of 
freight  to  cars,  the  composition  (in  terms  of  cars)  of  trains  leaving  a 
railyard,  train  routing,  the  scheduling  of  engines  to  train  routes,  the 
capacity  and  location  of  railyards,  and  many  other  aspects  of  rail  opera- 
tions. Modeling  is  complicated  by  "blocking"  or  "grouping"  contingencies 
in  which  cars  destined  for  several  final  locations  are  grouped  together  and 
treated  as  a unit  along  their  route  to  some  intermediate  railyard.  In  most 
instances,  the  number  of  possibilities  for  blocking  and  "reblocking"  is 
enormous . 

Although  analytic  approaches  have  been  successful  for  dealing  with 
certain  aspects  of  freight  flow  management  (Assad  (8]  delineates  efforts 
in  rail),  there  remains  great  potential  for  modeling  and  algorithmic 
development. 
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